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ABSTRACT 

A method of general perturbations utilizing Chebyshev series is used 
to investigate motion in the vicinity of the L 4 triangular point of the earth- 
moon system. The model used is that of the very restricted four body 
problem for the earth- moon-sun system. A harmonic orbit, in the numer- 
ical sense , with respect to a rotating coordinate frame centered at L 4 is 
found. The period of this harmonic orbit, being equal to the period of the 
disturbing force, is the same as the moon's synodic period. This orbit 
remains within 6860km of the L 4 point. It describes two different size 
loops about L 4 , the smaller one traversed in 36 percent of the period. The 
disturbing force, being nearly periodic with half the moon's synodic period, 
gives rise to another orbit about L 4 which is nearly periodic with half the 
synodic period of the moon. This orbit remains within 4574km of the L 4 
point for twelve periods investigated. Deviations from the mid orbit dur- 
ing this time is less than 381km. 
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A HARMONIC L 4 ORBIT 
FOR 

THE VERY RESTRICTED FOUR-BODY PROBLEM 
1. INTRODUCTION 

In some recent literature interest has been shown in the problem of the in- 
fluence of the sun on motion close to the libration points of the earth- moon sys- 
tem as well as motion about the earth and the moon itself. One possible model 
for the earth-moon-sun system in which the problem might be considered has 
been proposed by Su-Shu Huang (1960), who called it the "very restricted four- 
body problem. " Here the earth and moon describe circular orbits relative to 
one another, and their center of mass describes a circular orbit around the sun; 
all these orbits are Keplerian, lie in a plane, and no perturbations are con- 
sidered. Using this model Huang studied the motion of a fourth body of an in- 
finitesimal mass in a similar manner as in the restricted three-body problem. 
He concludes this model gives a general idea of where the fourth body could or 
could not go under given initial conditions when they are no longer very near the 
earth. Using a similar model Cronin et al. (1964) proved that under certain 
conditions the fourth body has a periodic motion, relative to a rotating coordi- 
nate frame, near each of the libration points of the restricted three-body prob- 
lem. Their proof is based upon assumptions concerning the masses and dis- 
tances of the bodies which are not satisfied by the earth-moon-sun system. 
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Siferd (1965) used Huang's model for the earth- moon-sun system to gener- 
ate some periodic orbits. Using a numerical integration procedure, the equa- 
tions of motion for the very restricted four-body problem were iterated upon 
utilizing a digital computer until some periodic orbits were obtained. By this 
technique eight periodic orbits, in the numerical sense, with a respect to a 
rotating coordinate system were found. Three periodic orbits were around 
the earth, three were around the moon, and two were around the earth- moon 
libration point (1^). No periodic orbits near the triangular points were 
obtained. 

Danby (1965) investigated the influence of the sun on motion close to the tri- 
angular points of the earth-moon system. He felt the very restricted four-body 
model inadequate for his investigation and therefore used a model in which the 
secular perturbations of the moon due to the sun were retained. The results 
may be said to strengthen the hope that stable motion around the triangular 
points of the earth-moon system is possible. Other investigators include Tapley, 
et al. (1963 and 1965) who used a model similar to the very restricted four-body 
model except the moon's orbit is inclined with respect to the earth-sun plane. 

The equations of motion for a particle near the triangular points of the earth- 
moon system are numerically integrated on a digital computer for various ini- 
tial conditions. One result indicates that a particle placed initially at a triangu- 
lar point (L 4 ) with zero relative velocity has an envelope of motion, centered at 
L 4 , going through a mode of expansion to a value of one earth-moon distance for 
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the major axis followed by a mode of contraction to a value of 1/8 earth-moon 
distance for the major axis. The envelope repeats this sequence several times 
during the 2500-day period investigated. The nature of these data suggests that 
such a motion may persist for a period of time much longer than that considered 
in the study. 

The present paper uses the very restricted four-body problem model as 
proposed by Huang for the earth-moon-sun system. The merits of this model 
for the earth-moon-sun system are still in doubt; however, it has been used in 
this study as a first attempt to find orbits which remain near triangular points 
of the earth-moon system. Using a technique proposed by Carpenter (1966), a 
harmonic orbit in the numerical sense , with respect to the rotating coordinate 
frame and about the L4 triangular point of the earth-moon system is found. In 
addition, a nearly periodic orbit having half the period of the harmonic orbit is 
obtained. 


2. THE MODEL AND COORDINATE SYSTEM 
Consider an infinitesimal body of mass, m in a system of three bodies m 1# 
m 2 , and m 3 which are the sun, earth, and moon respectively. Further assume 
that all four bodies remain in a plane so arranged that the center of mass, B of 
m 2 and m 3 is revolving around the center of mass, 0, of the entire system in a 
circular orbit and m 2 and m 3 themselves are revolving around B also in circular 
orbits. Table I indicates the numerical values used for this model, and Figure 
1 shows the geometry. 
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Physical Characteristics of the Model and System 
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The right angle X, Y axis system with its origin at the center of mass of the 
earth is rotating at a uniform rate so as to keep the center of mass of the moon 
on the X axis. The masses m, m! , m 2 , and m 3 are in the X, Y plane. A point 
L 4 , in the X, Y plane, is 60° from the X axis at lunar distance, a 3 , and in 
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advance of the moon’s position. This point corresponds to a triangular point of 
the earth- moon system three body problem. 

3. THE HARMONIC ORBIT 

Using the model described the motion of an infinitesimal mass, m in the 
vicinity of the L 4 point was investigated with the purpose of obtaining a harmonic 
orbit in the numerical sense. This harmonic orbit has a period equal to the 
period of the disturbing force, which for this model is the moon’s synodic period. 
Musen’s (1963) method is applied with the perturbations represented in Cheby- 
shev series as proposed by Carpenter (1966). In this method the geocentric 
position vector, "r, of the mass m near the L 4 triangular point is given by 

~r = (1 +ajr 0 +/3 w 

where a and /3 are the components of the perturbations, ~r 0 is the position vector 
in the fixed reference ellipse, a is the semimajor axis of the reference ellipse 
and 

-*• 1 dr 0 

w = TT 

n dt 

n being the mean motion in the reference ellipse and t the time. The functions 
a and /3 can be represented by uniformly convergent series in the interval -1 < 
r < 1 by 

Q ( T ) = a k T k (t) 
k=0 
00 / 

£(t) = X £k T k (t) 

k = 0 


6 



where the prime on the summation sign is used to indicate that the first term is 


to be factored by one-half. The T K (r) are the Chebyshev polynomials defined by 

T k (r) = cos[keos -1 r] 

where the coefficients of these polynomials are given by a k and /3 k . 

The synodic period for the model being utilized is given by the equation 



n 3 _n l 


where n, and n 3 are the mean motions of the sun and moon respectively. The 
dimensionless time r is related to time t from epoch by 

r = — 

P 

where t is zero at the epoch which is defined as the first time the mass n^ 
crosses the positive X axis. 

Starting with initial conditions at r = 0 in the X, Y plane and near the L 4 
point, an orbit was generated by using the Chebyshev polynomial procedure. 
Initial conditions corresponding to r of -1 and 1 were thus obtained. Using nu- 
merical partial derivatives from this orbit an iteration scheme was employed to 
match the initial conditions at r of -1 and 1. After several iterations this was 
accomplished, and the results are shown in Table n. The a' and j3' values, 
shown in this table, are derivatives of a and (3 with respect to nt. Relative geo- 
centric errors in position and velocity are indicated by the differences in Table 
II. These differences indicate agreement in ten significant figures which corre- 
spond to changes from r = -1 to r = 1 of 0.2 meters in the position vector, 
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TABLE II 

Initial Conditions at r = -1, 0, 1 for the Harmonic Orbit 


T 

a(r) x 10 3 

/3(t) x 10 3 

a'(r) x 10 3 

(3'( t ) x 10 2 

-1 

0 

1 

8.90721044 

4.41178027 

8.90721028 

1.43123468 

16.95842648 

1.43123517 

7.84005032 

14.7369698 

7.84005038 

-1.341882564 
-0.507181064 
-1. 341882534 

Differences 

a(+l) - a(-l) = - 1 . 6 x 10 -10 a’(+l) - a»(-l) - 6 . 0 xl 0 ' n 

/3(+l) -yQ(-l) = 4.9xlO ~ 10 £’(+1) -yS’(-l) = 3.0xl0 _1 ° 


rt and 5 x 10 " 7 meters per second in the velocity vector, m From the numer- 
point of view it seems adequate to call this a harmonic orbit. 

Chebyshev coefficients for this orbit are given in Table III. Using the initial 
conditions at epoch this orbit was extended for a total of twelve synodic periods 
(-1 < r < 25). This was done as a further check to insure the accuracy of the 
orbit. For this time period agreement with the initial orbit was eight significant 
figures in position and velocity. This is within the anticipated accuracy of the 
calculation. It was decided to plot the harmonic orbit with respect to a rotating 
coordinate system centered at L 4 . Referring to Figure 2, a geocentric vector, 
r 0 , directed toward L 4 has a magnitude, a, defined by 
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TABLE HI 


Chebyshev Coefficients for the Harmonic Orbit 


k 

a k x 10 6 

Pu * 10 6 

0 

9233.8688 

271.0928 

1 

-466.7498 

4715.9669 

2 

4193.3102 

-2590.7864 

3 

-2839.4208 

-2764.2264 

4 

1530.9949 

7775.0747 

5 

4927.4893 

-3235.6624 

6 

-1877.9548 

-5000.6993 

7 

-1940.7497 

1552.0981 

8 

507.1676 

1271.4117 

9 

356.8288 

-301.5204 

10 

-67.8275 

-172.1400 

11 

-41.2120 

36.8027 

12 

4.4834 

13.0026 

13 

4.3726 

-3.8665 

14 

0.1851 

-0.0154 

15 

-0.6623 

0.4549 

16 

-0.0905 

-0.1980 

17 

0.1213 

-0.0468 

18 

0.0042 

0.0448 

19 

-0.0196 

-0.0019 

20 

0.0048 

-0.0069 

21 

0.0024 

0.0029 

22 

-0.0021 

0.0007 

23 

-0.0001 

-0.0009 

24 

0.0005 

0.0000 

25 

-0.0001 

0.0002 

26 

-0.0001 

0.0000 

27 

0.0000 

0.0000 
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Figure 2. Coordinate Systems near L4 


LO 


where a is the semimajor axis of the reference orbit, kg is the earth’s gravita- 
tional parameter and n is equal to the mean motion of the moon, n 3 . The a, 

P coordinate system has its origin at r 0 with a directed along r^ and P at right 
angles to a in the direction of motion in the reference orbit. 

A S a, P coordinate system has its origin at (1 + a 0 ) r 0 which corresponds to 
L 4 . The value of a 0 is given by the quantity (a 3 - a)/a where a 3 , the moon's 
semi-major axis, is taken to be 60 earth radii. 

The 8a component is directed along r 0 and p is the same component previ- 
ously defined but translated parallel to itself to the L 4 point. Utilizing the data 
given in Table HI the harmonic L 4 orbit is plotted, see Figure 3, using the 8a, 

8a 
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P coordinate axes. Units for this plot are dimensionless with respect to the 
semimajor axis of the reference orbit. The harmonic orbit remains close to 
the L4 point having a deviation of less than 1.8 percent (6860 km) of the earth- 
moon distance. This orbit is seen to describe two different size loops about L 4 , 
the time in the smaller loop is 36 percent of the period. Conversion to the usual 
x, y coordinate system where x and y are centered at L 4 parallel to 8a, /3 but 
have dimensionless emits with respect to the semimajor axis of the moon is 
accomplished by the transformation 


and 


x = 8 a — — 
a 3 


Since the ratio (a 3 /a) is near unity there would be no noticeable difference by 
replacing 8a and (3 by x and y respectively in Figure 3, this matter being brought 
to attention for purposes of calculation. 


4. A NEARLY PERIODIC ORBIT 

Since the disturbing force is nearly periodic with half the synodic period of 
the moon, it is possible to find orbits which are nearly periodic with half the 
moon’s synodic period. One such orbit was obtained by approximately matching 
initial conditions at r of -1 and zero. The initial conditions are shown in Table 
IV. If the orbit were periodic with half the moon’s synodic period, initial con- 
ditions at t of -1 and zero would match with those at r = 1. The agreement in 
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TABLE IV 

Initial Conditions at r- -1, 0, 1 for the Nearly Periodic Orbit 


T 

a(r) x 10 3 

/S(r)xl0 3 

a’(r) x 10 2 

/3'( r ) x 10 3 

-l 

6.67701926 

9.31381378 

1.132682180 

-9.25579538 

0 

6.67701357 

9.31377754 

1.132682017 

-9.25577168 

l 

6.65506332 

9.24996973 

1.128066502 

-9.18511281 

Differences 

a(0) - a(-l) = -5.69xl(T 9 

a'(0) - a'(-l) = - 

■1. 63x 10 -9 

/3(0) -/3(-l) = -3.62x 10 -8 

/3'(0) -/3'(-l) - 

2. 37 x 10 -8 

a(0) - a(+l) = 2.19xl0" 5 

a’(0) - a'(+l) = 

4.615 x 10 -5 

/3(0) -/3(+l) = 6.38xl0~ 5 

/3'(0) ~{3'(+l) = - 

■7.06x 10~ 5 


initial conditions between r of -1 and zero is eight significant figures in position 
and velocity; however, between r of zero and one there is only five significant 
figure agreement. Further reduction of this difference between r of -1 and zero 
tended to increase the differences between zero and one . A plot of this nearly 
periodic L 4 orbit in the 8a, /3 coordinate system for r between minus one and 
one is given in Figure 4. Also shown is the envelope of the orbit for six synodic 
periods (-6 < r < 6). During this time interval the orbit remains within 1.2 
percent (4574 km) of the earth-moon distance from L 4 . Deviations from the r 
of minus one to one orbit are seen to be less than one tenth of a percent (381km) 
of the earth-moon distance. During longer time intervals the deviations are 
expected to increase. 


13 











8a 



Using the Chebyshev polynomial approach, motion in the vicinity of L 4 can 
be investigated for other models of the earth-moon-sun system. One of the 
more interesting models would include the moon moving in an eccentric orbit. 
With this model some insight may be gained as to the importance of the moon's 
eccentricity on motion near L 4 . The method is not restricted to simple models, 
e.g. , it is possible to study motion near L 4 using the actual motions of the prin- 
cipal bodies taken from a general theory or from an ephemeris. 
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Appendix A 


Derivations 


1. Derivation of system of equations for matching initial condition at r of -1 
and 1. 

Let a(-l), /3(-l), a' (~1), /3' (-1) and a(l), (6(1), a'(l^/3'(l) be the values 
at r = -1 and +1 respectively. If small changes are made in each of the pa- 
rameters at t = -1 these will result in changes in the parameters at r = +1 
and the ratios will give 

3a( l') da( 1) da( 1 ) da( 1 ) d/3( 1 ) d/3( 1 ) 

da(-l) ’ 3/3(~l) ’ da ( - 1) ’ d/3 (“1) ’ cia(-l) > 3/3(~l) > etc - 

for a total of 16 numerical partial derivatives. Now suppose a(-l> ft-1) , 

a'(-l) , and/^-l) are replaced by a(-l) + 8a , /3(— 1) + 8/3 , a'(-l) + Sa' , 

/3'(-l) + 8/3' then to first order the new values at t - 1 will be 

da(l) da( 1 ) 3a(l) Ba(l) 

a ( X ) + Ba(-l) 8a + 3/8(-l) + 3a'(-l) Sa> + 3/8' (-1) S Z 3 ' 




+ .Mu 

Ba(-l) 


Sa 


+ m i ). 

^(-i) 


S/3 + 


^(1) 

da' (— 1 ) 


Sa' 


, W) 


S/3' 


a'(l) + 


da ' ( 1 ) 
da(~l) 


Sa 


da ' ( 1 ) 

3/8(-l) 


S/3 + 


da ' ( 1 ) 
3a'(-l) 


Sa' 


3a' (1) 

W) 


8/3' 
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+ KH 

3a(~l) 


8a 


+ Am> s/3 + jgm 

3/3(~l) 3a' (- 1 ) 


Sa' 


, M(l) 


8/3' . 


The object is to find 8a, 8/3, 8a' , 8/3' such that the values at r = ±1 are equal. 
These 8 corrections will be the solution of the following system of equations 


1 3x(l) 

3a(~l) 


Sa - 


Ml) 

M~l) 


8/3- 


3a(l) 

3x'(-l) 


Sa' - 


Ml) 

3/3' (-1) 


8/3' 


a(l) - a(~l) = Ax 


Ml) 

3a(-l) 


Sa + 


Mi) 

M-i) 


8/3- 


Mi) 

3a'(-l) 


Sa' - 


Mi) 

M(-i) 


8/3' 


/3(1) -/3(-l) - A/3 


3a'(l) 3a'(l) 

6a 

3x(-l) 3/3(-l) 


3a'(-l) 


^'(1) 

M(-i) 


8/3' 


a'(l)-a'(-l) = Ax' 


3a(~l) 3/3(~l) ctx' (- 1 ) 


MOL" 

M(-i) 


8/3' 


/3'(1) - /S'(-l) - Zy3' 


In the above system of equations the unknown are 

8a , 8/3, Sa' , and 8/3' thus the equations were inverted numerically to find 
these unknows. Denoting the matrix of partial derivatives by P where 


3a( 1) 

3a(l) 

3a( 1 ) 

3a( 1) 

3a(~l) 

3/3(-l) 

3a ' ( - 1 ) 

M(-i) 

3/3(1) 

3/3(1) 

3/3(1) 

3/3( 1 ) 

3a(-l) 

3/3(-l) 

3a' (~1) 

M(-i) 

3a ' ( 1 ) 

3a'(l) 

3a ' ( 1 ) 

3a ' ( 1 ) 

3a(-l ) 

M-i) 

3a ' ( - 1 ) 

3/3'(-l) 

3/3'(l) 

3/3'(l) 

3/3' (1) 

M(i) 

3a(~l ) 

M-i) 

3a' (-1) 

3/3'(-l) 
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and define the matrices 



Sa 


Aa 


S/3 


A/3 

[SA] - 

Sa' 

[AA] - 

Aa' 


S/3 ' 


A/3' 


The system of equations can be written as 

(I - P) [SA] = [AA] 


where I is the unit matrix. 

The desired inverse of this system of equations is therefore 

[SA] = (I " P)" 1 [AA] 

For the harmonic orbit the P matrix was found to be 


- 6.519783 

-0.233560 

-1.091679 

- 3 

-19.707221 

0.9242190 

4.682355 

-10 

- 0.773883 

0.297094 

0.436769 

- 0 

12.07645 

0.37297 

1.72651 

6 


2. Derivation of the velocity equation 

The position vector, r , of a body is given by 


707652 

352294 

225023 

93592 


where a, /3, y are the components of the perturbations 7 0 , is the position 
vector in the fixed reference ellipse, a is the semimajor axis of the reference 
ellipse, R is the unit vector in the direction of the angular momentum of the 
motion in the reference ellipse and 


w 


i £0 

n dt 


n being the mean motion of the reference ellipse and t the time. Let prime 
denote derivative with respect to nt and take the derivative of the position 
vector, thus 


r' “ a'r 0 + /3‘ w + y' aR + (1 + a) 7' + /3w' 


Since w = 7' and w' 


1/n 2 d 2 7^t 2 where Keplerian motion is given by 


dt 2 



-n 2 a 3 



then w' - -(a/r 0 ) 3 7 0 making these substitutions 

7 0 + (/S' + 1 + a) w + y' aR 



3. Converting relative geocentric errors in position and velocity into errors in 
meters and meters per second. From Table II the errors were found to be 


Aa = 

- 1.6 x 

io - 10 

Aa' 

= 6.0 x 

10' 11 

A/3 = 

4.9 x 

10 -10 

A/3' 

= 3.0 x 

0 

1 

O 

tH 
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For motion in the plane the position vector is given by 


r - (l+a)7 0 +/3w 

in which 7 and w have the magnitude of the semimajor axis of the reference 
orbit, a = 59.756099826 earth radii. The position error vector is given by 

A? = Aa r 0 + A/3w 

which has a magnitude of 
| AT 1 = a j/ (Aa) 2 + (A/3) 2 

= (59.756099826)(6378165. ) /(-1.6 x 10" 10 ) 2 + (4.9 x 10 -10 ) 2 

= 0 . 196 meters 

For motion in the plane the velocity vector is given by 



7 0 + (/3' + 1 + a) w 


since the reference orbit is circular a = r Q and the velocity error vector is 
given by 

A v‘ - (A a' - A/3) r 0 + (A/3' + Aa) w 

which has a magnitude of 

|nAr' | = an / (Aa' - A/3) 2 + (A/3' + Aa) 2 


21 


From Table I 


n - 9.6595278 x 10 3 rad per hr. 

= 2.683202 x 10” 6 rad per sec. 

[(59. 756099826) (6378165)] 

[2.683202 x 10" 6 ] /[(0.6 - 4.9) x 10" 10 ] 2 + [(3.0 - 1.6) xlO' 10 ] 2 


4.62 x 10 -7 meters per second . 
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Tables of Chebyshev Coefficients 
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Table B-2 

Nearly Periodic Orbit for 6,26246 Synodic Periods (4400 hours). 


k 

a fc X 10 6 

X 10 6 

k 

a k X 10 6 

/ 3 k X ltf 

0 

6963.9090 

2606.5141 

1 

1267.7887 

- 873.0461 

2 

585.5660 

1302.5396 

3 

1363.4846 

- 963.0635 

4 

485.7971 

984.3795 

5 

1507.1638 

- 1119.4409 

6 

291.5963 

210.0762 

7 

1619.0140 

- 1318.8197 

8 

- 44.4945 

- 372.7762 

9 

1401.6762 

- 1331.1329 

10 

- 485.9954 

- 1189.4363 

11 

973.6196 

- 968.5940 

12 

- 880.8550 

- 1840.0052 

13 

204.6290 

- 116.2743 

lb 

- 927.9831 

- 2103.0428 

15 

- 844.5835 

860.5722 

16 

- 447.8515 

- 1367*0757 

17 

- 1659.5763 

1278.7061 

18 

3 l 4 . 93 ul 

788.3247 

19 

- 1141.2482 

895.9795 

20 

916.8582 

2348.2617 

21 

563.7213 

- 407.1751 

22 

620.6477 

1447.3271 

23 

1703.5065 

- 1421.5062 

24 

- 551.9072 

- 1326.9619 

25 

561.2978 

- 450.2972 

26 

- 951.8150 

- 2343.6191 

27 

- 1636.8243 

1334.8778 

28 

336.1607 

627.0519 

29 

- 800.7523 

654.1480 

30 

1014.6487 

2492.3961 

31 

1895.4409 

- 1547.6910 

32 

- 699.5445 

- 1714.7432 

33 

- 86.2864 

70.5438 

34 

- 616.3778 

- 1513.0319 

35 

- 1942.5146 

1588.8578 

36 

1369.8313 

3353.7686 

37 

2412.2744 

- 1974.8528 

38 

- 1241 . 49 ul 

- 3042.0659 

39 

- 1756.7531 

1435.0543 

40 

757.6500 

1861.6671 

41 

932.447 2 

- 761.0120 

42 

- 356.3634 

- 873.4677 

43 

- 390.6775 

321.0836 

44 

137.8553 

333.1282 

46 

136.0663 

- 113.2943 

46 

- 43 . 38 U 0 

- 105.6421 

47 

- 41.2262 

33.0096 

48 

10.1023 

28.8901 

49 

9.6678 

- 6.0441 

t >0 

- 1.8861 

- 7.9697 

51 

- 1.5876 

0.4239 

52 

0 • b 858 

1.4722 

53 

1.2066 

- 2.0122 

64 

- 1.4698 

0.7220 

55 

- 0.6766 

3.0527 

66 

2.4475 

0.0990 

57 

- 0.8579 

- 2.2271 

68 

- 2.4662 

- 0.9416 

59 

0.7778 

1.5086 

60 

1.1769 

0.0209 

61 

0.6425 

- 1.4536 

62 

0.0358 

1.0218 

63 

- 1.0284 

1.0925 

64 

- 0.0981 

- 0.6620 

65 

- 0.1118 

- 0.1593 

66 

- 0.4879 

- 0.6467 

67 

1.2122 

- 0.5766 

68 

0 . b 546 

1.0474 

69 

- 1.1221 

0.5551 

70 

- 0.1296 

- 0.2827 

71 

0.1560 

- 0.0122 

72 

- 0.6853 

- 0.8301 

73 

0.8686 

- 0.5366 

74 

0.9276 

1 .4402 

75 

- 1.4591 

0.8118 

7 b 

- 0.6324 

- 1.3886 

77 

1.4867 

- 0.7753 

78 

0.6022 

1.0401 

79 

- 1.0035 

0.5336 

80 

- 0.4319 

- 0.6962 

81 

0.4786 

- 0.2828 

82 

0.2578 

0.3899 

83 

- 0.2935 

0.1519 

84 

- 0.0663 

- 0.1624 

85 

0.1653 

- 0.0761 

86 

0.0776 

0.0779 

87 

0.0130 

0.0276 

88 

0.0055 

- 0.0249 

89 

0*0864 

- 0*0112 

90 

0*0243 

r 0.0070 

91 

0.0668 

. 0*0057 

92 

0 . 0 Z 03 

* 0*0043 

93 

0*0755 

“OT 0003 

94 

U . U 4 u 3 

- 0.0008 

95 

0.0544 

0*0079 

96 

0.0662 

- 0.0161 

97 

- 0.0038 

0*0026 

98 

0.0542 

- 0.0173 

99 

- 0.0954 

- 0.0182 

100 

0.0074 

- 0.0018 

101 

- 0.0958 

- 0.0132 

102 

- 0.0796 

0.0422 

103 

0.1021 

0.0162 

104 

- 0.0135 

0.0087 

105 

0.1468 

0.0236 

106 

0.1400 

- 0.0895 

107 

- 0.2996 

- 0.0586 

108 

- 0.1851 

0.1075 

109 

0.2457 

0.0530 

110 

0.0866 

- 0.0719 

111 

- 0.1260 

- 0.0321 

112 

- 0.0489 

0.0329 

113 

0.0475 

0.0206 

114 

0.0108 

- 0.0111 

115 

- 0.0126 

^• 0*0077 

116 

0.0209 

0.0028 

117 

0.0032 

0.0065 

118 

0.0141 

- 0.0007 

119 

- 0.0008 

- 0.0067 

120 

0.0063 

- 0.0001 

121 

- 0.0009 

- 0.0013 

122 

- 0.0197 

0.0001 

123 

- 0.0004 

- 0.0017 

124 

- 0.0190 

0.0002 

125 

0.0006 

0 . 0046 

126 

0.0032 

- 0.0001 

127 

0.0000 

0.0033 

128 

0.0251 

- 0.0000 

129 

- 0.0002 

- 0.0028 

130 

- 0.0068 

- 0.0000 

131 

- 0.0006 

- 0.0053 

132 

- 0.0250 

0.0002 

133 

0.0005 

0.0025 

134 

- 0.0098 

0.0000 

135 

0.0008 

0.0060 

136 

0.0283 

- 0.0002 

137 

- 0.0006 

- 0.0042 

138 

0.0014 

- 0.0001 

139 

- 0.0007 

- 0.0051 

140 

- 0.0315 

0.0002 

141 

0.0007 

0.0081 

142 

0.0223 

- 0.0001 

143 

- 0.0001 

- 0.0014 

144 

0.0125 

-o.oooi 

145 

- 0.0005 

- 0.0076 

146 

- 0.0390 

0.0003 

147 

0.0012 

0.0120 

148 

0.0437 

- 0.0003 

149 

-o.ooii 

- 0*0111 

160 

- 0.0338 

0.0002 

151 

0.0004 

0.0077 

162 

0.0207 

- 0.0001 

153 

- 0*0003 

- 0.0043 

154 

- 0.0103 

o.oooi 

155 

0.0005 

0-0020 

156 

0.0044 

- 0.0001 

157 

-o.oooi 

- 0*0009 

158 

- 0.0018 

- 0.0000 

159 

- 0.0003 

0.0003 

160 

0.0006 

0.0001 

161 

0.0000 

- 0.0001 
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Table B-3 


Harmonic Orbit for Seven Synodic Periods. 


k 

<x k X 10* 

,ys k X 10® 

k 

a k X 10® 

/3 k X 10® 

0 

7457,2159 

810.6086 

1 

- 146.2737 

1921.1129 

2 

1205.6185 

- 500.8064 

3 

- 224.7309 

1623.8637 

4 

1394.2414 

- 722.6099 

5 

- 428.5244 

924.6166 

6 

1664.0341 

- 761.5695 

7 

- 760.0051 

- 245.9927 

8 

1463.3128 

- 48.6037 

9 

- 966.9615 

- 1597.7711 

10 

841.1339 

1826.0628 

11 

- 394.3421 

- 2190.0430 

12 

- 138.0494 

3797.1057 

13 

1405 . 5604 

- 1068.2085 

14 

- 573.9662 

2569.3025 

15 

2931.6451 

587.8964 

16 

46.4183 

- 3005.8293 

17 

1049.7570 

- 477.3923 

16 

- 166.0784 

- 4544.2480 

19 

- 2081.8566 

- 1803.0256 

20 

- 1769.4628 

2083.0325 

21 

119.1789 

3110.2725 

22 

824.9464 

- 2196.4307 

23 

- 2058.7499 

34.2711 

24 

661.8066 

3171.5108 

25 

678.7966 

32.8783 

26 

704.3463 

1176.6787 

27 

1565.6424 

- 1448.8345 

28 

- 672.0417 

- 1272.6793 

29 

442.0562 

- 330.1177 

30 

- 667.4497 

- 2159.4794 

31 

- 1649.9648 

1347.2723 

32 

460.5920 

1134.1503 

33 

- 479.6221 

398.4775 

34 

673.7834 

2146.8321 

35 

1872.3926 

- 1531.3793 

36 

- 840.0863 

- 2060.5426 

37 

- 515.0913 

420.1151 

38 

- 340.4363 

- 634.8463 

39 

- 1636.3100 

1254.2673 

40 

1228.2742 

3014.9956 

41 

2346.5570 

- 1923.5473 

42 

- 1285.7716 

- 3160.8588 

43 

- 1930.7972 

1586.8700 

44 

682.4341 

2168.2507 

45 

1140.5773 

- 932.6490 

46 

- 461.4496 

- 1130.6934 

47 

- 533.3931 

432.2833 

48 

194.4513 

477.1465 

49 

206.5843 

- 171.3310 

50 

- 68.3608 

- 170.3621 

51 

- 66.9316 

60.4684 

62 

20.9647 

32.6676 

53 

19.4166 

- 15.7319 

64 

- 6.0197 

- 13.2520 

55 

- 5.8660 

- 1.1356 

56 

2.3297 

2.8936 

57 

0.7112 

3.6051 

68 

- 1.1327 

- 1.5764 

59 

0.6108 

1.2709 

60 

- 0.7226 

0.3364 

61 

0.6086 

- 6.6775 

62 

2.0755 

0.8289 

63 

- 0.7801 

0.0019 

64 

- 1.5674 

- 0.0627 

65 

- 0.5831 

- 5.7733 

66 

0.4467 

- 0.9704 

67 

0.9449 

3.3120 

68 

- 0.3348 

0.4266 

69 

0.2556 

- 2.1731 

70 

0.7828 

0.7733 

71 

- 1.0809 

1.5228 

72 

- 0.6363 

- 0.8429 

73 

0.4054 

- 0.5742 

74 

- 0.1065 

- 0.2232 

75 

0.8127 

- 0.3245 

76 

0.6327 

1.0732 

77 

“ 1.1764 

0.5902 

78 

- 0.4940 

- 0.8443 

79 

0.4232 

- 0.2123 

80 

- 0.0828 

- 0.1386 

81 

0.6938 

- 0.3732 

82 

0.6253 

1.0494 

83 

- 1.3982 

0.7573 

84 

- 0.8525 

- 1.4063 

85 

1.4817 

- 0.0174 

86 

0.7809 

1.2646 

87 

- 1.1785 

0.6598 

88 

- 0.6671 

- 0.9125 

89 

0.7856 

- 0.4349 

90 

0.3412 

0.5694 

91 

- 0.4601 

0.2382 

92 

- 0 * 17 o 6 

- 0.3133 

93 

0.2342 

- 0.1113 

94 

O .0776 

0.1444 

95 

- 0.0942 

. 0.0542 

96 

- 0*0438 

“ 0*0467 

97 

■'010224 

- 0*0346 

98 

0.0291 

0.0049 

99 

o.oooi 

0*0192 

100 

— O.Olul 

0.0003 

101 

0.0044 

- 0.0005 

102 

- 0.0075 

0.0099 

103 

- 0.0158 

- 0.0108 

104 

O .0126 

- 0.0198 

105 

0.0227 

0*0108 

106 

- 0.0088 

0.0230 

107 

- 0.0228 

- 0.0073 

108 

0.0066 

- 0.0205 

109 

0.0183 

0.0075 

110 

- 0.0092 

0.0151 

111 

- 0.0124 

- 0.0098 

112 

0.0109 

- 0.0096 

113 

0.0073 

0.0088 

114 

- 0.0072 

0.0053 

115 

- 0.0039 

- 0.0034 

116 

-0.0000 

- 0.0027 

117 

0.0019 

- 0.0026 

118 

0.0053 

0.0012 

119 

- 0.0009 

0.0051 

120 

- 0.0053 

- 0.0005 

121 

0.0003 

- 0*0030 

122 

0.0012 

0.0001 

123 

- 0.0000 

- 0.0012 

124 

0.0038 

- 0.0001 

125 

0.0000 

0.0049 

126 

- 0.0071 

0.0002 

127 

- 0.0002 

- 0.0067 

128 

0.0080 

- 0.0002 

129 

0.0002 

0.0066 

130 

- 0.0071 

0.0002 

131 

- 0.0000 

- 0.0054 

132 

0.0053 

- 0.0000 

133 

- 0.0001 

0.0038 

134 

- 0.0035 

- 0.0001 

135 

0.0001 

- 0.0023 

136 

0.0020 

0.0000 

137 

0.0000 

0.0013 

138 

- 0.0010 

0.0001 

139 

- 0.0001 

- 0.0006 

140 

O . 00 o 5 

- 0.0001 

141 

0.0001 

0*0002 

142 

- 0 . 00 O 2 

0.0001 

143 

- 0.0000 

- 0.0001 

144 

O.OOol 

0.0000 

145 

- 0.0001 

0-0001 

146 

- O.OOol 

- 0.0001 

147 

0.0001 

- 0.0001 

148 

O.UOOl 

0.0001 

149 

-o.oooi 

0.0001 

150 

- O.OOul 

- 0.0001 

151 

0.0001 

- 0.0001 

152 

O.OOol 

0.0001 

153 

- 0.0001 

0.0001 

154 

- 0 . 00 O 0 

- 0.0000 

155 

0.0001 

- 0.0000 

156 

0 . 00 O 0 

G .0000 

157 

-o.oooi 

0.0000 

158 

- 0.0000 

- 0.0000 

159 

0.0000 

0.0000 
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Appendix C 

Computer Programs 


Table C-l 

Listing of CHBY, Subroutines, and a Sample Case. 

SJOB 1203P003 126 M32RK M3081C1234 
{EXECUTE IB JOB 

SIBJOB SOURCE. MAP. GO 

SIBFTC CHBY DECK.M94.XR7 

MODIFIED FOR THE FOUR BODY PROBLEM 
***«*• MOON MUST BE THE FIRST DISTURBING BODY CONSIDERED IN THIS 
VERSION OF THE PROGRAM ****** 

CARPENTER. CHBY PROGRAM. COMPUTES PLANETARY PERTURBATIONS (ALPHA* 
BETA* GAMMA) IN CHEBYSHEV SERIES 

DOUBLE PRECISION EKX,XJ»V.SA»SB»SE»SN.SGZ»SM.EP»P.Q.R»EIC.RAD» 

IP I *TWOPI»EP2»S02.SI2»C02»SGZZ.SMINNR.SO»sI.CO»SN2»A.B»C»THET.DLL* 
2E.COSC.SINC.ESE.ECE,BK.DCOS»DSIN.COS2 C»SiN2C,X»VAL.VALZ»FAC»ZK* 
3TWOSE.CN1»CN2*CN3»CN4.DN1»DN2»CNN,SAB»TRaTIO»DXZR»DSQRT*TEPOCH. 
4PEPOCH.DTCHEB.PTCHEB.EPSl 

COMMON X»EkK.NTRM.Nl.XJ.V.SA.SB.SE.SN»SGZ«SM*EP*P»Q»R.IPRINT*NTIM 
COMMON TRAT IO.DXZR.ZX 

DIMENSION X< 701.8) *AA<6)»ZM6>* V(5 ) • I TITLE! 12) 

DIMENSION A(3) »B( 3) *C< 3 ) *P( 3 ) *Q( 3 ) *R< 3 ) *vAL( 3 .2 *3 ) »VALZ( 3*2 ) 

1 FORMAT (4F15.7) 

2 FORMAT (1H08X1HX17X4HTHET41X1HV/1P7D18.10) 

3 FORMAT ( 80( 1H3 ) ) 

4 FORMAT (34H0 THE MOTION OF G IN DTCHEB DAYS IS. 

411PD25.15.8H DEGREES) 

5 FORMAT (I8.6P6F12.4) 

6 FORMAT (24H0VALUE OF SERIES AT X"-l /1P6 d20*12/ 

61 24H0VALUE OF SERIES AT X> 0 /1P6D20.12/ 

62 24H0VALUE OF SERIES AT X**l /1P6D20.12) 

7 FORMAT ( 26H0 VALUE OF SERIES AT EPOCH/ UX1HA19X1HB19X 
71 1HC16X7HA PRIME13X7HB PRIME13X7HC PRIME//) 

8 FORMAT ( 1P6D20. 12 ) 

9 FORMAT (/2(7X1HX1X11HALPHA*1.E6 2XlOHBETA*l*E6 1X11HGAMMA*1.E6 )// 
91) 

10 FORMAT <2<I8<6P3F 12*4)) 

12 FORMAT ( 10H CARPENTER) 

13 FORMAT (22H0 INTEGRATION CONSTANTS/ 1H0.9X2HZK. 1 1 »5< 17X2HZK* 1 1 > > 

13 FORMAT (918) 

17 FORMAT ( 8H0 TEPOCH12X6HDTCHEB12X6HPEPOCH12X6HPTCHEB13X4HEPS1/5D18. 
Ill) 

22 FORMAT I 17H0TEST ZERO VALUES/1P6D20.12) 

23 FORMAT ( 5H0 X55X2HN113X2HN213X2HN312X4HCOSE1 1X4HSINE// ) 

26 FORMAT < I3.48X.1P3E13.7) 

29 FORMAT (26H0 INITIAL CONDITIONS AT X» .F3.0/11X1HA19X1HB19X 
291 1HC16X7HA PRIME13X7HB PRIME13X7HC PRIME//) 

32 FORMAT < 7H0VAL ID .F10.1.6H THRU .F10*l) 

33 FORMAT (4D18.11) 

34 FORMAT < I 8.0P3F 12.0 . 10X . 1P3E13.7) 

36 FORMAT (48HJ N1 NTRM ITRM IPRlNT IPUNCH IC) 

37 FORMAT <2 8H0REFERENCE ELLIPTIC ELEMENTS/ 

371 1H013X2HSA23X2HSE22X3HS0222X3HSI2/1PAD23.15) 

38 FORMAT < 1HO12X3HCO223X2HSM22X4HSGZ220X6HSMINNR/ 1P4D23* 13 ) 

39 FORMAT < 1H013X2HSB22X3HSN2/1P2D25. 15 ) 

41 FORMAT < 12A6) 

43 FORMAT <1H1) 

52 FORMAT (64H0VECTORIAL CONSTANTS REFERRED TO EQUATOR AND MEAN EQUIN 
5210X 1950.0/ 

522 1H010X1HP19XIHQ19XIHR19X1HA19XIHb19X1HC/< 1P6D20. 12) I 

PI*3. 141592653589793 
TWOPI=2.*PI 
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Table C-l (Continued) 


RAD=TW0PI/360. 

EK IS THE GAUSSIAN CONSTANT 
EK=0. 017202098950+00 

THE NEXT EK. IS FOR EARTH! REF • ITEM) UNITS FOR EK SQ ARE ER CUBE/HR SQ. 

E K=4. 46 19969 18D+0Q 

EP2 IS NEWCOMB OBLIQUITY FOR JAN 0 1950* 23DEG 26MIN 44.84SEC 

EP2=23. 44578888888889 
EP=EP2*RAD 

ENTRY POINT FOR SUCCESSIVE CASES 

55 CONTINUE 
DO 59 1=6*8 
DO 59 J=l» 701 
59 X <J.I)=0. 

READ (5,41) ITITLE 

ITITLE IS 72 BCD CHARACTERS USED AS A TITLE FOR THE RUN 
READ (5,15) N1,NTRM,ITRM,IPRINT*IPUNCH,IC 
Nl+1 NUMBER OF SPECIAL VALUES 

NTRM NUMBER OF COEFFICIENTS (NOT TO EXCEED 701) 

I TRM NUMBER OF ITERATIONS FOR THIS RUN 

IPRINT =1 PRINT RESULTS ONLY 

=2 PRINT SPECIAL VALUES AND EXPANSION COEFFICIENTS 
=3 PRINT TEST QUANTITIES IN FFD SUBROUTINE 
I PUNCH =1 DO NOT PUNCH COEFFICIENTS ON CARDS 
=2 PUNCH COEFFICIENTS 
IC =1 INITIAL CONDITIONS AT X=-l 

=2 INITIAL CONDITIONS AT X= 0 

=3 INITIAL CONDITIONS AT X=+l 

READ (5,33) T EPOCH *DTCHEB * EPS1 

TEPOCH JULIAN DATE OF THE EPOCH 

DTCHEB TIME INTERVAL IN DAYS FOR VALIDITY OF THEORY 

EPS1 CRITERION FOR TRUNCATING COMPUTED SERIES (AFTER STATEMENT 410) 

READ (5,33) P EPOCH *PTCHEB 

PEPOCH, PTCHEB ARE EPOCH AND TIME INTERVAL FOR DISTURBING PLANETS 

ASSIGN 570 TO NPUNCH 
IF ( IPUNCH-U62. 62*60 
6u PUNCH 12 
PUNCH 3 
PUNCH 3 
PUNCH 3 

IF (ITRM-1) 61,61,62 

61 ASSIGN 568 TO NPUNCH 

62 CONTINUE 

READ (5,33) SA.SE ,S02 »SI 2 »C02 »SM» SGZ2 ,SMl NNR 
C 

C NOTATION FOR ELEMENTS OF DISTURBED PLANET 
C SA SEMI-MAJOR AXIS IN EARTH RADII 

C SE ECCENTRICITY 
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Table C-l (Continued) 

C S02 ARGUMENT OF PERIHELION (DEGREES) 

C SI INCLINATION (DEGREES) WITH RESPECT TO ECLIPTIC 

C A VALUE OF 23.445788888 PUTS BODY IN ECLIPTIC PLANE 

C C02 LONGITUDE OF ASCENDING NODE (DEGREES) 

C SM MASS OF DISTURBED PLANET IN UNITS OF EARTH MASS 

C SGZ2 MEAN ANOMALY AT EPOCH (DEGREES) 

C THE NEXT CARD IS NOT NEEDED WHEN EARTH IS PRIMARY 
C SMINNR MASS OF INNER PLANETS (TO BE ADDED TO SOLAR MASS) 

C 

READ (5,33)ZK. 

ZK ARE INTEGRATION CONSTANTS USED FOR COMPUTING SPECIAL VALUES 

S0=S02*RAD 
SI =S I 2*RAD 
C0=C02*RAD 
SGZ=SGZ2*RAD 

SN = EK*USQRT ( 1 . +SMI NNR+SM ) /SA**1 • 5 
EKX=.5*SN*DTCHEB 
TRAT 10=DTCHEB/PTCHEB 
DXZR = < T EPOCH-P EPOCH ) *2 . /PTCHES 
SN2=SN/RAD 

SB=SA*DSQRT( l.-SE**2) 

SAB=SA/SB 
WRITE (6.43) 

WRITE ( 6 *41 ) I T I TLE 
WRITE (6,36) 

WRITE (6,15) N 1 »NTRM , I TRM ,IPRINT»I PUNCH , IC 
WRITE (6,37) SA ,SE ,S02 , S 1 2 
WRITE (6,38) C02 , SM , SGZ2 .SMINNR 
WRITE (6,39) SB »SN2 

WRITE (6,17 ) T EPOCH , DTCHEB , P EPOCH, PTCHEB , EPS1 
X J1 = TEPOCH-DT CHEB/2 . 

XJ2=TEPOCH+DTCHEB/2. 

WRITE (6,32) XJ1.XJ2 
BK=2.*EKK/KAD 
WRI TE ( 6 ,4 )BK. 

VALZ ARE THE INITIAL CONDITIONS 

READ (5,33) ( < VALZ < I , J ) , I =1 , 3 ) , J= 1 ,2 ) 

X J= I C-2 

WRITE (6,29) XJ 

WRITE (6,8) ( ( VALZ ( I , J ) ,1=1, 3), J= 1.2) 

WRITE (6,13) ( I , 1=1 ,6) 

WRITE (6,8) ZK. 

PQRDP2 

CALL PURDP2 (SO, SI ,CO,P,Q,R,EP) 

XJ = 0. 

CALL PQRDP2(SO,SI , CO ,P , Q ,R , XJ ) 

DO 65 1=1,3 
A ( 1 ) = SA*P ( I ) 

B ( I )=SB*Q( I ) 

65 C(I)=SA*R(I) 

WRITE (6,52) <P< I ) ,Q( I ) ,R< I ) ,A( I ) ,B( I ) ,C( I ) ,1=1 ,3) 

READ INPUT PERTURBATIONS 

NT I M=0 
WRITE (6,9) 

70 READ (5,5) K • ( AA ( I ) , I = 1 ,6 ) 
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Table C-l (Continued) 


EACH CARO CONTAINS THE COEFFICIENTS OF T SUB I OF X AND T SUB K+l 
OF X FROM THE PREVIOUS ESTIMATES OF ALPHA* BETA* GAMMA 

KP=K+1 

WRITE 16.10)K.(AA< I ) » I«1 .3 ) *KP» < AAI1) *1-4*6) 

IF (K-1001) 80*92*92 
80 K1=K+1 
KP=K+2 

IF (Kl-700) 82*82*70 
82 CONTINUE 
NTIM*KP 
DO 90 1*1*3 
J*I+3 

X ( K1 ♦ 1+5 ) *AA < I ) 

X ( XP * I +5 ) *AA ( J ) 

90 CONTINUE 
GO TO 70 
92 CONTINUE 
I TYPE*1 

CALL FFD TO READ DISTURBING PLANET DATA 

CALL FFD (ITYPE) 

1TTRM*1 

ENTRY POINT FOR SUCCESSIVE ITERATIONS 

94 ITYPE*2 

DO 96 1*1*5 
DO 96 J*l. 701 
96 X (J.I)*0. 

THE FOLLOWING OPERATIONS THRU STATEMENT 230 61 VE THE CHEBYSHEV 
EXPANSIONS OF THE DISTURBING FORCES. THE EXPANDED QUANTITIES ARE 
V(l) V ( 2 ) VI 3 1 VC4) VI 5 ) 

M1*ADR 2»*M2*ADR M3*AOR COSCE) SIN(E) 

WHERE Ml. M2. M3 ARE THE SUBSCRIPTED QUANTITIES FROM THE EQUATIONS 
FOR THE PERTURBATIONS. ADR IS THE FACTOR SMALL A/ SMALL R SUB ZERO* 
AND E IS THE ECCENTRIC ANOMALY IN THE REFERENCE ELLIPSE 


XJ=-l.D+0 

CALL FFD FOR FIRST SET OF SPECIAL VALUES 

CALL FFD (ITYPE) 

THET=PI 

IF (IPRINT-2U10. 100*100 
100 CONTINUE 

WRITE ( 6*2 )XJ .THET »<V( I ) *1*1*5) 

110 CONTINUE 

DO 200 1=1*5 
V< 1)=.5D+Q*V( I ) 

200 CONTINUE 

COEFF ADDS CONTRIBUTIONS OF SPECIAL VALUES TO SERIES COEFFICIENTS 


CALL COEFF (XJ.V) 

EN1=N1 

DLL=P I / EN 1 
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Table C-l (Continued) 

DO 210 X*2»N1 
THET-THET-DLL 
XJ*DCOS( THET ) 

C 

C CALL FFD FOR INTERIOR POINTS 
C 

CALL FFO ( ITYPE I 
C 

CALL COEFF IXJ.V) 

IF <IPRINT-2)210»202»202 
202 WRITE (6»2)XJ .THET . (V( I > .1*1.5) 

210 CONTINUE 
XJ*+l.D+0 
C 

C CALL FFO FOR LAST SET OF SPECIAL VALUES 
C 

CALL FFO ( ITYPE) 

IF ( I PR 1 NT-2) 2 14 *2 12 *2 12 
212 CONTINUE 
THET*0. 

WRITE (6.2 )XJ .THET » (V( I ) .1*1.5) 

214 CONTINUE 

00 220 1*1.5 
V( I )*.5D+0*V( I > 

220 CONTINUE 
C 

CALL COEFF (XJ.V) 

FAC*2./EN1 
N2*NTRM— 10 
00 230 X=1»NTRM 
00 230 1*1.5 
X(X.I)*X(X.I)«FAC 
230 CONTINUE 

IF ( I PRINT~2)364. 240.240 
240 CONTINUE 

WRITE (6*43) 

WRITE (6.41) ITITLE 
WRITE (6.25) 

DO 340 X=1 .NTRM 
XP*X-1 

WRITE (6.26) XP.(X(X,1) .1*1.5) 

340 CONTINUE 
C 

C TEST ZERO VALUES OF X(X,1) THRU X(K.5) 

C 

00 362 1*1.5 
ZM I )*.5*X(1,I ) 

BK*1. 

00 362 k*3.N2.2 

ZX ( I ) =ZM I ) +BX*X ( X . I ) 

362 CONTINUE 

WRITE (6.22) (ZX( I ) .1*1.5) 

C 

C PERFORM SERIES INTEGRATIONS (THRU STATEMENT 410). THE PREVIOUS 
C VALUES OF THE PERTURBATION COEFFICIENTS ARE NO LONGER NEEDED. THE NEW 
C SERIES OCCUPY THE OLD LOCATIONS (X(X.6)» X(K.7). X(X.8)) 

C 

364 TW0SE*2.*SE 

X(1»4)*X(1 .4) — TWOSE 
C 
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Table C-l (Continued) 


C MLTPLY (11.12.13) CAUSES THE SERIES X(X»U) AND X(X,I2) TO BE 
C MULTIPLIED AND THE PRODUCT TO BE STORED IN X(X.I3) 

C 

CALL MLTPLY (4.1.8) 

X( 1»4)=X( 1.4I+TW0SE 
CALL MLTPLY (5.2.7) 

DO 370 X= 1 *N2 
X(X.8)*X(X»8) “X ( X *7 )' 

370 CONTINUE 
C 

C NTGRT (11.12) CAUSES THE SERIES X(X.U) TO BE INTEGRATED WITH RESPECT 
C TO NT AND THE INTEGRAL TO BE STORED IN XIX. J2) 

C 

CALL NTGRT 18.7) 

CALL MLTPLY' (5.7.6) 

CALL MLTPLY (5.1.8) 

CALL MLTPLY (4.2.7) 

DO 380 X=l»N2 
X(X.8)=-X(X.8)-X(X.7) 

380 CONTINUE 

CALL NTGRT (8.7) 

X(1»4)*X( 1 .4 )— TWOSE 
CALL MLTPLY (4.7.8) 

X ( 1 .4 ) =X ( 1 .4 ) +TWOSE 
DO 390 X=1»N2 
X(X»6)=X(X»6)+X(X*8) 

X(X»1 )=-SE*X(X»4) 

390 CONTINUE 

X(l.l)=2.+X(l.l) 

CALL MLTPLY (1.2.7) 

CALL NTGRT (7.8) 

DO 400 X s 1 .N2 
X(X.6)=X(X.6)+X(X,8) 

X(X,8)=.5*X(X.8)-2.*X(X»6) 

400 CONTINUE 
C 

C ALPHA IS NOW IN X(X.6) 

CALL NTGRT (8.7) 

C 

C BETA IS NOW IN X(X.7) 

X ( 1 .4 ) ®X ( 1 .4) —TWOSE 
CALL MLTPLY (3.4.1) 

CALL NTGRT (1.2) 

CALL MLTPLY (2.5.8) 

CALL MLTPLY (3.5.1) 

CALL NTGRT (1.2) 

CALL MLTPLY (2.4.1) 

DO 410 K=1 »N2 
X(X»8) S X(X»8) “X ( X* 1 ) 

410 CONTINUE 
C 

C GAMMA IS NOW IN X(X.8> 

C 

C ESTABLISH TRUNCATION POINT. NTIM* FROM MAGNITUDES OF COEFFICIENTS 
C 

EPTST=EPS1/1.E1Q 
NT IM=NTRM+1 
DO 414 X= 1 .NTRM 
NT IM=NT IM-1 
DO 412 1=1.3 

IF ( ABS ( X (NTIM.I+5) )-EPTST) 412.416.416 
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Table C-l (Continued) 


412 CONTINUE 
414 CONTINUE 
416 CONTINUE. 

WRITE 16.43) 

WRITE ( 6.41 ) I T I TLE 

PRINT COEFFICIENTS 

418 CONTINUE 

WRITE (6.9) 

DO 460 X=1»NTIM.2 

X1=X-1 

L=X+1 

DO 430 1=1.3 
J=I+3 

AA( I ) =X (X » 1+5 ) 

AA( J )=X (L . 1+5 ) 

430 CONTINUE 

WRITE (6.101X1. (AA( I ) .1=1.3) .X.IAAII) .1=4*6) 

460 CONTINUE 

EVALUATE SERIES AT EPOCH 

DO 500 1=1.3 
DO 500 J= 1 .2 
DO 500 X=1 * 3 
VAL ( I .J.X)»0. 

500 CONTINUE 

DO 520 1=1.3 
IS* I +5 

VALl I .1.1 )«.5*X(1 . IS) 

VAL ( I. 1 .2 ) = • 5*X ( 1 . IS ) 

VAL( 1.1.3)=. 5*X(1. IS) 

VALl 1.2*1 )*0. 

VAL (1.2.2) =0. 

VAL (1.2.3) =0. 

BX=1. 

FX=1 • . 

DO 510 X=3.N2.2 
BX=-BX 
. FX= FX+2. 

VAL (1.1.1) *VAL ( 1.1.1) “X ( X— 1 . I S ) +X ( X . I S ) 

VAL( I .1.2 )*VAL( 1.1.2) +BX*X(X.1S) 

VAL( I .1.3)=VAL( 1.1 » 3 ) +X ( X-l *1S)+X(X»IS) 
CNl=(FX-2. )**2 
CN2=(FX~1.)**2 
CN3=— BX#( FX-2. ) 

VAL ( I »2»1)=VAL( 1 . 2 » 1 ) +CN1*X ( X-l .IS )—CN2*X ( X. IS ) 
VAL ( 1.2.2 )=VAL( 1.2 .2 ) +CN3*X ( X-l * IS ) 

VALI1.2.3) =VAL (1.2.3) +CN1*X (X-1.IS)+CN2*X(X.1S) 
510 CONTINUE 

VAL (I.2.1)=VAL(I.2.1) /EXX 
VALl I .2.2) =VAL ( I. 2. 2) /EXX 
VALl I .2.3>*VAL( 1.2. 3) /EXX 
520 CONTINUE 

WRITE (6.41) 1 TITLE 
WRITE (6.7) 

WRITE (6.6) (( (VALl 1 *J*X) .1 = 1.3) »J=1»2) »X=1» 3) 

X J= IC— 2 

WRITE (6.29) XJ 

WRITE (6. 8) ( (VALZI I *J) .1=1.3) .J=1.2) 
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Table C-l (Continued) 


CONSTANTS OF INTEGRATION FOR INITIAL CONDITIONS 
XJ*IC~2 

CN1«SGZ+EKK*XJ 

E*CN1 

S1NC*DSIN(E) 

COSC-DCOSIE) 

00 542 1*1*3 

ESE*SE*SINC 

ECE*1.-SE*C0SC 

BK.* ( CN 1 -E+ESE I / ECE 

E»E+BK-.5*BK*BK*ESE/ECE 

SINC*OSIN(E) 

COSC*OCOS(E) 

542 CONTINUE 

CN1*C0SC/U.-SE*C0SC) 

CN3*SINC/ ( 1 •-SE*COSC ) 

CN2—SINC 

CN4*C0SC-SE 

DN1*VALZ (3*1) — VAL ( 3* 1 • I C ) 

DN2*VALZ( 3,2 )-VAL( 3.2.10 
ZM5)=UN1*CN1+DN2*CN2 
ZM 6 ) *DN1*CN3+DN2*CN4 
DN1*VALZ(1.1)-VAL(1.1.IC) 

DN2*VALZ< 1, 2» -VAL( 1.2.10 

ZM 3 ) =2 .*DN1+VALZ < 2 . 2 ) -VAL ( 2 *2 . IC) 

DN1=DN1-2.*ZM3) 

ZM 1 )*UN1*CN1+DN2*CN2 
IK ( 2 )*DN1*CN3+DN2*CN4 
S I N2C=2 . *S 1 NC*COSC 
C0S2C»1.-2.*S1NC**2 
CNN*-2.+SE*SE 
DN1*VALZ(2.1)-VAL(2»1.IC) 

DN2*-3.*(SE*ZM 1KZM3) )*EKK*XJ 

ZM4)*DN1+DN2 +(-CNN*SINC-.5*SE*SIN20*ZM 1) 

1 ♦ < -2.*COSC+.5*SE*COS2C)*ZM2) 

WRITE ( 6* 13 ) ( 1 * 1*1 *6 ) 

WRITE (6* 8) (ZM I ) *1*1.6) 

ADO CONTRIBUTIONS FROM CONSTANTS TO SERIES COEFFICIENTS 

X(1*7)-X(1,7)+2.*ZM4) 

X(2*7)»X<2»7)+3.*CSE*ZM1)-ZM3))*EMC 

XU,6)*X<1.6)+4.«ZM3)-TW0SE*ZM1) 

X ( 1 *8 ) *X( 1 *8 ) -TWOSE*ZM 5 ) 

X( 1 .4 ) *X ( 1 ,4) +TWOSE 
CALL MLTPLY (4*4*1) 

X(1*1)*X(1*1)-1. 

X(K.l) IS «5*C0S( 2E) 

CALL MLTPLY (4.5,2) 

X(K,2) IS .5*SIN(2E) 

DO 550 K*1,N2 

X(M6)*X<M6)+ZM 1)*X(K,4)+ZM2)*X(K»5> 

X(M8)=X(M8) +ZM 5 ) *X(K»4 )+ZM6 )*X (K.5 ) 
X(M7)=X<K.7)+ZM1)*<CNN*X(K»5)+SE*X<M2) ) 

1 +ZM2)*( 2«*X(K,4)— SE*X(M1) ) 

550 CONTINUE 

PRINT AND PUNCH COEFFICIENTS 
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Table C-l (Continued) 


WRITE (6*43) 

WRITE (6*41 ) IT ITLE 
WRITE (6*9) 

DO 570 X«1*NTIM*2 

Xl-X-1 

L*X+1 

DO 564 1*1*3 
J«l+3 

AA(I)*X(X*I+5) 

AA< J)-X(L.l+5) 

564 CONTINUE 

WRITE (6»10)X1»(AA(I) »I*1»3) *X»(AA(I) .1=4*6) 

GO TO NPUNCH .(568*570) 

568 PUNCH 5 »X1 * (AA( I > .11*6) 

570 CONTINUE 

EVALUATE SERIES AT EPOCH (SHOULD EQUAL INITIAL CONDITIONS) 

DO 580 1*1 »3 
DO 580 J*l»2 
DO 580 X*1.3 
VAL(I*J»X)»Q. 

580 CONTINUE 

DO 600 1*1*3 
IS*I+5 

VAL(I,1,1)*.5*X(1,IS) 

VAL( 1.1*2 ) **5*X ( 1 * IS) 

VAL ( 1*1 *3 )=«5*X( 1 » IS ) 

VAL (1*2.1) *0. 

VAL( 1*2*2) *0* 

VAL(I.2»3)“0. 

BX*1. 

FX*1. 

DO 590 X*3.N2.2 

BX—BX 

FX* FX+2. 

VAL(I*1.1)*VAL( 1.1.1 )-X(X-l.IS)+X( X, IS) 

VAL( 1*1*2 )*VAL( 1*1*2) +BX*X(X,IS) 

VAL ( I • 1 » 3 ) * VAL ( I » 1 » 3 ) +X ( X-l * I S ) +X ( X ♦ I S ) 

CNl»(FX-2.)**2 
CN2*(FX-1. )**2 
CN3*-BX«(FX-2.) 

VAL( I * 2 *1 )*VAL( I *2*1 )+CNl#X( X— 1 * IS )-CN2*X (X* IS) 

VAL ( 1 .2 »2 ) »VAL ( I *2*2) +CN3*X { X-l* IS) 

VAL(I*2.3)*VAL(I *2.3)+CNl*X(X-l.lSI+CN2»X(X»ISI 
590 CONTINUE 

VAL( 1*2*1 )*VAL( 1*2*1) /EXX 
V AL ( I » 2 .2 ) « V AL ( I » 2 ♦ 2 ) / EXX 
VAL( I »2*3)*VAL( I *2 *3) /EXX 
600 CONTINUE 

WRITE (6*41) ITITLE 
WRITE (6*7) 

WRITE (6*6) (((VAL(I.J*X).I = 1*3).J*1.2)*X=>1.3) 

DO 602 1*1*6 
ZX(I)*0. 

602 CONTINUE 

ITTRM*1TTRM+1 

IF I I TTRM-I TRM) 94*610*604 
604 PUNCH 3 
PUNCH 3 
PUNCH 3 
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GO TO 55 

610 IF ( I PUNCH- 1 ) 94.94,620 
620 ASSIGN 568 TO NPUNCH 
GO TO 94 
END 

SIBFTC FFO DECK. *M94 *XR7 

SUBROUTINE FFD ( I TYPE ) 

C SUBROUTINE FOR COMPUTING SPECIAL VALUES OF DISTURBING FORCE 

DOUBLE PRECISION EKK ,XJ » V *SA .SB .SE »SN ,SGZ « SM, EP ,P ,U,R , EK.RAD »X « 
1SAB,P7,P9»P11»P13»SA2,CU,C31 * SAP » SEP ,SOP2 *SI P2 ,COP2 *SMP2 »S6ZP2 » 
2SMINNR,SOP,COP,SlP.SGZP,SNP,SNP2,PPP,QPP,RPP,SBP,AP.BP,CP,PP,QP, 
3RP.SMP»SMPA2 .EKP,ENT,SG.E»SNE»CSE,DCOS,DsIN.DSQRT , ESE , ECE .BK .COSE . 
4SINE.SR»SR2»SR3.ADR.TADR.SRZ.Sl»S2»EX»SRv»F»ENTP.SGP»CKl,CK2.CK3, 
5SPZ,RH02,SRP2,SPV,RH0,CRH0,CSRP.TRATI0 »Dx2R,XJP,ZK,CS1,CS2,SS2. 
6HALFSE »SEE »BARY »BR 
DOUBLE PRECISION PK1 , PK2 ,PK3 , SRZDS 

COMMON x,erk,ntrm.ni,xj»v»sa,sb.se,sn,sgz»sm,ep»p,q.r,iprint.ntim 
COMMON TRAT IO.DXZR »ZK 

DIMENSION SEP ( 9 J »SGZP (9) »PPP < 3) »UPP(3) »RPP(3) . AP ( 3 .9 ) * BP ( 3 » 9 ) * 
1CP(3*9) »PP ( 3 ) *OP ( 3 ) * RP ( 3 ) *SMPA2 ( 9 ) *X ( 701.8). Y(3,1200), AA(6). 

2NTPM19 ) *SRZ(2 ) »SW ( 2 ) ,S1 < 2 ) .S2 ( 2 ) .PERT ( 3 ) ,S<3) .SRVl 3) .F (3) .EKP(9) . 
3SPZ ( 3 ) >SWZ ( 3 ) .SPV ( 3 ) . RHO ( 3 ) . I T I TLE ( 12 ) .P ( 3 ) .Q ( 3 ) . R ( 3 ) ♦ V ( 5 ) . 
4NFIRSTI9) *2M6) .BR(3) 

1 FORMAT ( 12A6 ) 

2 FORMAT ( 1H1 . 12A6.4H M= . 1 2 ) 

3 FORMAT I4D18.il I 

4 FORMAT ( 12HJSMALL 0ME6AF 15 . 7 .F13. 7 . 15H ) ECLIPTIC AND5X7HSMALL E 

41 F12.8/12H CAP OMEGAF 15. 7 . F 13. 7 . 15H ) MEAN EQU1NOX5X7HSMALL A 

42 F15.7.6X2HER/8H SMALL IF 19 .7 .F 13 . 7 . 11H > 1950.09X7HSMALL N 

43 E16.8.8H DEG/HR /7H G ZER0F20. 7 , FI 3 . 7 »20X7HSIflALL ME16.8) 


5 FORMAT 

(1HK9X1HA16X1HB16X1HC13X6H P 1QX6H 

Q 10X7H 

R // ) 

6 FORMAT 

(6F17.8.15H 

) 

EQUATOR AND/ 



61 

6F17.8.15H 

) 

MEAN EQUINOX/ 



62 

6F17.8.12H 

) 

1950.0) 



7 FORMAT 

(1H0) 





8 FORMAT 

(6F17.8.15H 

) 

P.Q.R. SYSTEM/ 



81 

6F17.8, 9H 

) 

OF/ 



82 

6F17.8.15H 

) 

COORDINATES) 



9 FORMAT 

( I8.6P6F12.4) 




10 FORMAT 

(/2(7X1HK1X11HALPHA*1.E6 2X10HBETA*1»E6 

IX 1 1HGAMMA*! • E6 ) 6X 


1017HST0RAGE// ) 

11 FORMAT (58H0 K A**2*F (I ) SUM ALPHA, BETA, GAMMA E.COSIEl.SIN 

111(E) 6X6HSRV ( I ) 12X6HSRZ ( I ) 13X5HSW ( I ) ) 

12 FORMAT (1HO»I3»1PD16.10,E18.7.5D18.10> 

13 FORMAT (4X.1PD18.10.E18.7.5D18.10) 

14 FORMAT (25H0 M DISTURBING PLANE TS94X6HRHO (I) ) 

16 FORMAT (2( I8.6P3F12.4) ,218) 

17 FORMAT ( 1H017X7HDEGREES7X7HRADI ANS ) 

18 FORMAT <6H SNP =1PD24. 15 .6HRAD/HR ) 

IF ( ITYPE-2 ) 100.300.300 

COMPUTE CONSTANTS AND READ DISTURBING PLANET DATA 

100 CONTINUE 
KEP I TR=3 

IF (SE-0. 25)104. 102. 102 
102 KEP I TR*4 
104 CONTINUE 

RAD=2. *3. 141592653589793/360. 

SAB=SA/SB 

C THE NEXT EK IS FOR EARTH(REF. ITEM) UNITS FOR EK SQ ARE ER CUBE/HR SQ. 


36 



n n n n n n n n 


Table C-l (Continued) 


EK=4 #4619969180+00 
P 7 = 7./ 6. 

P 9= 9./ 8. 

Pll*ll./10. 

P13*13./12. 

SA2*SA*SA 

HALFSE=SE*.5 

SEE=-2.+SE**2 

C11*1.5*SA2 

C31=-15.*SA2/8. 

M*Q 

K2=G 

ENTRY FOR SUCCESSIVE DISTURBING PLANETS 

110 CONTINUE 
M=M+1 

RE AO <5.1)1 TITLE 
WRITE(6.2)ITITLE»M 

****** MOON MUST BE THE FIRST DISTURBING BODY CONSIDERED IN THIS 
VERSION OF THE PROGRAM ****** 

READ <5.3) SAP .SEP (M ) *S0P2 .SIP2.COP2 .SMP2 .SGZP2 .SMINNR 

ELEMENTS OF DISTURBING PLANETS IN SAME FORM AS DISTURBED PLANET 

S0P=S0P2 *RAD 
C0P*C0P2 *RAD 
SIP=SIP2 *RAD 
SGZP<M)=SGZP2 *RAD 

SNP=ER*DSURT< 1.+SMINNK+SMP2 ) /SAP**1.5 
SNP2 =SNP /RAD 

WRITE (6,17) 

WRITE <6. 4) SOP 2 .SOP .SEP < M) . C0P2 .COP .SAP ,Sl P2 ,S I P .SNP2 . SGZP2 . 
1SGZP <M ) .SMP2 
WRITE (6,18) SNP 

CALL PURDP2 < SOP .SI P .COP ,PPP .QPP . RPP » EP ) 

SBP=SAP *DSQRT < l.-StP<M)**2) 

DO 120 1*1.3 

AP < I »M ) =SAP *PPP<I) 

BP<I,M)*SbP *QPP ( I ) 

CP < I »M ) *SAP *RPP ( 1 ) 

120 CONTINUE 
WRI TE ( 6 .3 ) 

WRITE(6.6) (AP( I »M) »BP< I »M) »CP< I ,M) ,PPP ( I ) ,QPP( I ) ,RPP< I ), 1*1,3) 
DO 130 1*1.3 
PPtI ) *0 • 

QP < I ) =0. 

RP ( I ) =0. 

130 CONTINUE 

DO 140 J=l»3 

PP<1 )*PP<1 )+P(J)*PPP< J) 

QP < 1 ) *QP <1 ) +P < J ) *QPP l J ) 

RP < 1 ) =RP ( 1 )+P< J)*KPP(J) 

PP ( 2 ) *PP < 2 )+Q< J)*PPP< J) 

QP < 2 ) *QP < 2 )+Q< J ) *QPP ( J ) 

RP < 2 )*RP<2 )+Q( J)*RPP( J) 

PP ( 3 ) =PP ( 3 )+R< J)*PPP< J) 

QP < 3 ) *QP <3 ) +R < J ) *QPP < J ) 

RP < 3 ) =RP < 3 )+R ( J ) *RPP ( J ) 

140 CONTINUE 


37 


nnorinn nnn 


Table C-l (Continued) 


00 150 I = 1 » 3 

AP ( I »M ) =SAP *PP(I> 

bP l I »M ) =SBP *QP(I) 

CP ( I »M J =SAP *RP(I) 

150 CONTINUE: 

WRITEI6.7) 

WRITEI6.8 ) (AP( I »M) »BP( I *M) *CP< I »M) »PP( I ),QP(I ),RP(I ),I=1.3> 

SMP=SMP2/( l.+SMINNR+SMI 

SMPA2 ( M ) =SA2*SMP 

EK.P ( M ) =EK1C*SNP / ( SN*TRAT 10) 

REAO PERTURBATIONS. ALL COEFFICIENTS STORED IN Y ARRAY 

NFIRST <M) =K.2+1 
NTPM(M)*0 
WR I TE ( 6. 10 ) 

190 READ (5*9) 1C » < AA (I) * I =1 *6 ) 
tCK=K.+l 

KP=k+NFIRST(M) 

K.P1=ICP + 1 

WRITE (6. 16 )K, (AA( I ) » 1 = 1.3) *K.K»(AA( I > ♦ 1=4,6 ) .KP.KP1 
IF ( K— 1001 >200*110*220 
200 K.1=KP 
K2=K.P1 

IF (K.2-1200I202. 202.190 
20 2 CONTINUE 

NTPM(M)=Kw+2 
DO 210 1=1.3 
J= I +3 

Y(I.ICl) =AA( I ) 

YII.K.2) *AA ( J ) 

210 CONTINUE 
GO TO 190 
220 NPLNET=M 

BARY=SMPA2 ( 1 ) / ( SA2+SMPA2 ( 1 ) ) 

K .=0 

RETURN 

COMPUTE SPECIAL VALUES 

300 CONTINUE 
tC=K+l 


COMPUTE POSITION VECTOR OF DISTURBED PLANET 

ENT =EK.K.*X J 
SG=SG2+ENT 
E=SG 

SNE=DSIN(SG ) 

CSE=DCUS(SG ) 

DO 310 1=1 .K.EPITR 

ESE=SNE*SE 

ECE=1.-CSE*SE 

BK.= ( SG -E+tSE)/ECE 

E=E+BK-.5*BK.*BK*ESE/ECE 

CSE=DCOS(E) 

SNE=DSIN( E ) 

310 CONTINUE 
CQSE=CSE 
SINE=SNE 
CS1=C0SE-SE 


38 



nnn r\nn n r> n 


Table C-l (Continued) 


CS2-2.*C0SE-SE*< .5-SINE**2 ) 

SS2-SEE*SINE+SE*(SINE*C0SE+3.«ENT> 

V<4)«C0SE 

VI5I-SINE 

SR -SA*<1.-SE*C0SEJ 

SR2«SR*SR 

SR3-SR*SR2 

AOR-SA/SR 

TA0R»2.*ADR 

SRZ(l)-SA*CSl 

SRZi2)*SB»SlNE 

SW C 1 ) » -S A*ADR* S 1 N E 

SW<2>« SB*ADR*COSE 

S1U)- AUR*COSE 

Sl<2)« AOR*SlNE*SAb 

$2(1)= -TAOR*SINE 

S2 ( 2 ) - TA0R*CS1»SAB . , 

ALPHA-2 .*ZK ( 3 > +ZK. < 1 ) *CS1 +Z.K. ft ) »S I NE 
GAMMA- ZM5)*CS1+ZM6)*SINE 

BETA »ZIC ( 4 ) +ZM 1 ) *SS2+ZK. C 2 ) *CS2-3 . *ENT*ZK ( 3 ) 

NT IM IS THE NUMBER OF TERMS IN THE SERIES FOR THE PERTURBATIONS 

EX»2.*XJ 

IF (NTIM-1J370. 330*340 
330 ALPHA-. 5*X< 1*6) +ALPHA 

BETA *« 5*X (1*7) +BETA 

GAMMA*. 5*X< 1*8) +GAMMA 
GO TO 370 
340 CONTINUE 

00 360 IS-6 *8 
B2-0. 

Bl-0. 

NTR-NTIM-l 
DO 350 I-l.NTR 
Il-NTR-1+2 

BZ=EX*bl-B2+X< II. IS) 

B2*ol 

Bl-BZ 

350 CONTINUE 
JS-IS-5 

PERT ( JS )».5*(EX*B1-2.*B2+X 1 1 * IS ) ) 

360 CONTINUE 

ALPHA-PERT ( 1 ) +ALPHA 

BETA -PERT < 2 ) +BETA 

GAMMA-PERT < 3 ) +GAMMA 

370 CONTINUE 

S < 1 ) «ALPHA*SRZ < 1 > +bt TA*SW ( 1 ) 

S < 2 ) -ALPHA*SRZ < 2 ) +BE TA*SW < 2 ) 

S ( 3 ) -GAMMA* SA 

EXPAND SOLAR TERMS IN POWERS OF PERTURBATIONS 

SQ=S ( 1 ) *S ( 1 ) +S < 2 ) *S ( 2 ) +S ( 3 ) #S ( 3 ) 

SRV(1)=SRZ<1)+S<1> 

SRV ( 2 ) *SRZ ( 2 ) +5 ( 2 ) 

SRV( 3 ) - S<3) 

DEL-<2.*ISRZ(1)*S(1)+SRZ<2)*S<2) )+SQ)/SR2 
C1-C11*SQ/(SR2*SR3) 

C2-C11*DEL/SR3 

C3-C31*0EL**2*( l.-P7*OEL*( l.-P9*0EL*( l.-pll*DEL*( l.-P13*DEL) ) ) )/SR 
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13 

F ( 1 )=C1*SRZ( 1 )+C2*S( 1 )+C3*SRV ( 1 ) 

F ( 2 ) =C1*SRZ(2 )+C2*S< 2 )+C3*SRV ( 2 ) 

F ( 3 ) = C2*S(3)+C3*SRV(3) 

********** F(1),F(2)*F(3) CHANGED TO THOSE SHOWN BELOW******* 
PK1 »PK2 »PK3 »SRZDS » F 1 * F2 *F3 ARE DOUBLE PRECISION 
PK1=SA2/(SR3*( l.+DEL)**1.5) 

PK2=SA2/SR3 

PK3=3.*PK2/SR2 

SRZDS=SRZ( 1 )*S( 1)+SRZ(2 )*S(2 > 

F<1 )=(-PKl+PK2)*SRV( 1)-PK3*SRZ< 1)*SRZDS 
F ( 2 ) = ( — PK1+PK2 ) *SRV (2)-PK3*SRZ(2) *SRZDS 
F ( 3 ) = ( -PK1+PK2 I *SR V ( 3 ) 


371 


372 


IF { IPRINT-3) 372,371*371 
WRITE (6,11) 

WRITE (6,12) K»F ( 1 ) * ALPHA * E»SRV(1) 
WRITE (6,13) F ( 2 ) , BETA ,COSE ,SRV ( 2 ) 
WRITE (6,13) F<3) , GAMMA, SINE, SRV ( 3 ) 
CONTINUE 


,SRZ( 1) ,SW ( 1 ) 
»SRZ ( 2 ) »SW ( 2 ) 
»SRz ( 3 ) »SW ( 3 ) 


NEXT TAKE THE DISTURBING PLANETS 


IF < IPRINT-31376,374,374 
374 WRITE (6,14) 

376 XJP=XU*TRAT IO+DXZR 
BR( 1 ) = 0« 

BR ( 2 ) -0 • 

BR ( 3 ) =0 • 

EX=2 ,*X JP 

DO 480 M=1 ,NPLNET 

ENTP=EKP(M)*XJP 

SGP=SGZP ( M ) +ENTP 

E=SGP 

SNE=DSI N( SGP ) 

CSE=DCOS< SGP ) 

DO 380 1=1,3 
ESE=SNE*SEP(M) 

ECt=l.-CSE*SEP(M) 

BK= ( SGP-E+ESE ) /ECE 
E=E+BK-.5*BK*BK*ESE/ECE 
CSE=DCOS( E ) 

SNE=DS I N( E ) 

380 CONTINUE 

CKl=CSt-SEP(M) 
CK2=-SNE/(1.-SEP(M)*CSE) 
CK3=CSt/ ( 1 •— SEP (M ) *CSE ) 

DO 390 1=1,3 

SPZ( I )=AP( I ,M ) *CK1+8P ( I ,M)*SNE 
SWZ( I ) =AP ( I ,M)*CK2+BP( 1 ,M)*CK3 
390 CONTINUE 

IF (NTPM(M)-l 1400,410,420 
400 ALP=0. 

BEP=0. 

GAP=0. 

GO TO 450 
410 JPP=NF I RST ( M) 

ALP= ,5*Y ( 1 , JPP ) 

BEP= ,5*Y ( 2 , JPP ) 

GAP=,5*Y ( 3, JPP ) 

GO TO 450 
420 CONTINUE 
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NTP=NTPM(M)-1 
DO 440 IS*1 *3 
B2=0. 

bi=o* 

DO 430 1=1 *NTP 
1 1=NTP“ 1+2 
JPP=NF I RST (M)+I 1-1 
BZ=EX*bl-B2+Y( IS.JPP) 

B2=B1 

bi*bz 

430 CONTINUE 

JPP=NF IRST (M) 

PERT ( IS )=*5*( EX+B1-2. *B2+Y ( I S» JPP ) ) 

440 CONTINUE 

ALP=PERT ( 1 ) 

BEP=PERT ( 2 ) 

GAP=PERT < 3 ) 

450 CONTINUE 
RH02=Q • 

SRP2=0 • 

DO 460 1=1*3 

SPV( I )=(1.+ALP)*SPZ( I )+BEP*SWZ(I)+GAP*CP( I»M)+BR< I ) 
RHO(l)=SPV( I >-SRV< I > 

RH02=RH02+RH0 ( I ) *RHO { I ) 

SRP2=SRP2+SPV( I )*SPV( I ) 

460 CONTINUE 

CRH0=SMPA2 ( M ) / ( RH02+DSQRT ( RH02 ) ) 

CSRP=SMPA2 IM) / ( SRP2+DSQRT (SRP2) ) 

DO 470 1=1*3 
BR( I )=BARY*SPV< 1 ) 

Fill =F t I ) +CRHO*RHO( I ) -CSRP*SPV( I ) 

470 CONTINUE 

IF < IPRINT-3)480»471»471 

471 WRITE 16*12) M*F ( 1 ) »ALP * E.SPVU ) »SPZ< 1 ) »SWZ < 1 ) *RHO( 1 ) 

WRITE (6*13) F(2).BEP.CSE,SPV<2)»SPZ<2)»SWZ<2)*RHO(2) 

WRITE 16.13) F(3) *0AP»SNE,SPV(3) *SPZ(3) »SWZ<3) »RH0(3) 

480 CONTINUE 

V(l) =F(1)*S1< 1)+F(2)*S1(2) 

V ( 2 ) =F( 1)*S2( 1)+F(2)*S2(2) 

V ( 3 ) =F( 3 ) 

RETURN 

END 

JIBFTC PQRDP2 NODECK. »M94»XR7 

SUBROUTINE PQRDP2 ( SO »SI .CO.P .Q.R.EP ) 

C COMPUTE VECTORS P,Q*R 

DOUBLE PRECISION SO . Si » CO.P *0 »R »EP . EPC .EPS.CCO »SCO »CSO *SSQ. 
X CSI »SSI .TP1.TP2.DC0S.DSIN 
DIMENSION P ( 3 ) »Q(3) »R(3) 

C EP IS OBLIOU1TY OF ECLIPTIC 
EPC=DCOS(EP) 

EPS=DSIN(EP) 

CCO=DCOS( CO) 

SCO=DS I N( CO ) 

CSO=DCOS( SO ) 

SSO=DS I N ( SO ) 

CSI =DCOS (SI) 

SSI =DS 1 N ( S I ) 

P ( 1 ) =-CSI*SSO*SCO+CSO*CCO 
TP1 =+CS I *SSO*CCO+CSO*SCO 
TP2 =+SSI*SSO 
P(2)=EPC*TP1-EPS*TP2 
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P(3)=EP$*TP1+EPC*TP2 

Q ( 1 ) *-CS I *cso*sco-sso*cco 

TP1 = CSI*CSO*CCO-SSO*SCO 

TP2 = SSI*CSO 

Q(2)=EPC*TP1-EPS*TP2 

Q(3)=EPS*TPi+EPC*TP2 

R(l> = SSI*SCO 

TP1 =-S5I*CC0 

TP2 = CSI 

R(2)=EPC*TP1-EPS*TP2 

R(3)=EPS*TP1+EPC*TP2 

RETURN 

END 

SIBFTC COEFF DECK..M94 *XR7 
SUBROUTINE COEFF 

C COMPUTE CONTRIBUTIONS OF SPECIAL VALUES TO SERIES COEFFICIENTS 
DOUBLE PRECISION XJ.V.X. TZ.T1.T2 *TwOX*EkR 
DIMENSION V ( 5 J *X( 701.6) 

COMMON X *EXK.*NTRM*N1 »XJ » V 
TZ = 1 • 

T 1 = XJ 

DO 100 1=1*3 
X IlflUX ( 1 * I ) +V ( I ) 

x ( 2 *n=x ( 2 *n+v( i )*ti 

100 CONTINUE 
TW0X=2.*XJ 
DO 120 K= 3 .NTRM 
T2*TW0X*T 1-TZ 
DO 110 1=1.3 
X (K.n=X U»I)+V(2>*T2 
110 CONTINUE 
TZ=T 1 
T1 = T2 

120 CONTINUE 
RETURN 
END 

SIBFTC MLTPLV DECK.M94,XR7 

SUBROUTINE MLTPLY (11*12*131 
C PERFORM CHEBYSHEV SERIES MULTIPLICATION 

C THE FACTORS ARE IN X(K*J1J AND XOC.J2J. THE PRODUCT GOES IN XU.J3). 
COMMON X( 701*8 ) *EKK»N 
DOUBLE PRECISION EKK.X.EX 
J1 = U 
J2= I 2 
J 3* I 3 

DO 100 K= 1 * 701 
X <K * J3 ) *0* 

100 CONTINUE 
N 1 = N 

DO 130 K*1*N1 
EX =0. 

IM*N-K.+ 1 

DO 120 1 = 1 * IM 

L=I+K-1 

EX » EX -t-X( I.J1)*X(L.J2)+X(L.J1)*X( I*U2) 

120 CONTINUE 

X ( K. * J3 ) *• 3*EX 
130 CONTINUE 

XU*J31=XU.J3)-.5*X(1,J1)*X(1*J2) 

DO 150 (C=3.N1 
EX s 0* 

I M=K-1 

DO 140 I =2 » IM 
LaK-i+1 

EX * EX +X( I. J1)*X(L.J2| 

140 CONTINUE 

xtx.j3 )=.5 *ex ♦x(K.J3) 

150 CONTINUE 
RETURN 
END 

SIBFTC NTGRT DECX.M94 *XR7 

SUBROUTINE NTGRT (J1.J2) 

C PERFORM CHEBYSHEV SERIES INTEGRATION 

C THE INTEGRAND IS STORED IN X(K.Jl) AND THE INTEGRAL IN XIK*J2) 

C THE FACTOR EKK COMES FROM D( NT ) * EICK * D(X) 

COMMON *( 701.8) *EKK*N 
DOUBLE PRECISION EKK»X*EK2.EK 
I=J1 
J«J2 

X(1.J)«0. 

EK2=2./EXK 
EX = 0. 

DO 100 K = 2 * N 
EK=EK + EK2 

X(X.U)= ( X(K-*1 * I )— X( K+l * I ) ) / EK 
100 CONTINUE 

X { N+l * J ) = X ( N • I ) / ( EK.+EK2 ) 

RETURN 

END 
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Table C-2 

Listing of INVT, Subroutines, and a Sample Case. 

SJOB 1203P002 126 M32RK M3Q81C1234 
♦EXECUTE IB JOB 

S I B JOB SOURCE .MAP .GO 

SIBFTC INVT DECK..M94.XR7 

C 

DIMENSION CM (4*1) * BM (4*4) 

1 FORMAT (4018.1) 

23 FORMAT ( 1HL13X2HZ1 14X2HZ2 14X2HZ314X2HZ414X2HZ5 14X2HZ6/4X6P6F 16.8 > 

29 FORMAT ( 17H0N0RMAL EOUAT IONS/ IX ) 

30 FORMAT! 1P4E15. 7. 5X.1PE15.7) 

A3 FORMAT (1H1) 

41 FORMAT (12A6) 

READ (3.1) ( (BM( I»J),I=1»4)»J=1»4) 

READ (3*1) (CM(I»1)»I=1»4) 

WRITE ( 6 * 43 ) 

WRITE (6.41) ITITLE 
WRITE (6.29) 

WRITE (6. 30) ( (BM( I .J) »J=1 .4) ,CM( I ,1 ) .1=1,4) 

CALL MAT INV(BM»4»CM» 1 .DETERM ) 

Z1=CM( 1.1) 

Z2=CM( 2 * 1 ) 

Z3=CM (3.1) 

Z4=CM (4,1) 

WRITE (6,23) Z1»Z2»Z3,Z4,Z5,Z6 

RETURN 

END 

SIbFTC MAT INV DECX.M94.XR7 

SUBROUTINE MAT INV ( A.N ,B »M .DETERM ) 

DIMENSION IP I VOT (4) »A(4,4) » B ( 4 » 1 ) .INDEX (4*2) .PIVOT (4) 

INITIALIZATION 

10 DETERM=1.0 
15 DO 20 J=1.N 
20 I PIVOT ( J ) =0 
30 DO 550 1=1, N 

SEARCH FOR PIVOT ELEMENT 

40 AMAX=0.0 
45 DO 105 J=1.N 
50 IF ( IPIVOT(J)-l) 60.105,60 
60 DO 100 (C=1.N 

70 IF ( IPIVOT(K)-l) 80,100.740 
80 IF (ABS (AMAX)-ABS ( A ( J ,K ) ) ) 85.100,100 
85 IROW=J 
90 I C0LUM=K 
95 AMAX = A ( J « K ) 

100 CONTINUE , 

105 CONTINUE 

110 IP 1 VOT ( ICOLUM ) = I P I VOT ( I COLUM ) +1 

INTERCHANGE ROWS TO PUT PIVOT ELEMENT ON DIAGONAL 

130 IF ( IROW-ICOLUM) 140,260,140 
140 DETERM=-DETERM 
150 DO 200 L= 1 »N 


48 



nn n nnn n n n 


Table C-2 (Continued) 

160 SWAP=A< IROW.L) 

170 A ( I ROW »L ) =A ( I COLUM. L ) 

200 A! ICOLUM, L)=SWAP 
205 IFIM) 260.260.210 
210 00 250 L=1.M 
220 SWAP=B ( I ROW »L ) 

230 B ( I ROW » L ) =B ( 1 C0LUM • L ) 

250 B( ICOLUM,L)*SWAP 
2 60 INDEX! 1.1 )*IR0W 
270 INDEX! 1.2 )*ICOLUM 
310 PIVOT! I )*A< ICOLUM. ICOLUM) 

DIVIDE PIVOT ROW BY PIVOT ELEMENT 

330 All COLOM . I COLUM ) - 1 . 0 
340 DO 350 L=1 »N 

350 A ( 1 COLUM , L ) *A ( I COLUM. L ) /P I VOT ( 1 ) 

355 IF ! M> 380.380.360 
360 DO 370 L= 1 »M 

370 B< 1 COLUM, L)*B( I COLUM. L) /PIVOT ( I ) 

REDUCE NON-PIVOT ROWS 

380 DO 550 L1=1,N 

390 IF! LI— ICOLUM) 400.550.400 

400 T = A ( L 1 » I COLUM ) 

420 A ( LI . I COLUM ) *0 • 0 
430 DO 450 L=1.N 

450 A!L1»L)*A(L1,L)-A( ICOLUM, L)*T 
455 IFIM) 550.550,460 
460 DO 500 L=1«M 

500 b!Ll»L)=B(Ll,L)-B( ICOLUM,L)*T 
550 CONTINUE 

INTERCHANGE COLUMNS 

600 DO 710 1=1, N 
610 L*N+1-I 

620 IF ( INDEX ( L , 1 ) — INDEX I L . 2 ) ) 630.710,630 
630 JROW* I NDEX ( L, 1 ) 

640 JCOLUM* I NDEX ( L .2 ) 

650 DO 705 K*1»N 
660 SWAP*A!K.»JROW) 

670 A ( X « JROW ) *A( K» JCOLUM ) 

700 A I K .JCOLUM ) 'SWAP 
705 CONTINUE 
710 CONTINUE 

DO 800 1*1 »N 
J=N+1-I 

800 DETERM*DETERM*P I VOT 1 J ) 

740 RETURN 
END 
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Table C-2 (Continued) 


S.DATA 


7.6519783 

00 

19.707211 

DO 

0.773883 

DO -12.07645 

DO 

A 


0.23356 

00 

0.073761 

DO 

-0.297094 

DO -0.37297 

DO 

b 


1.091679 

00 

-4.682333 

DO 

0.563231 

DO -1.72651 

DO 

A 

PRIME 

3.707632 

00 

10.332294 

DO 

0.225023 

DO -5.93592 

DO 

6 

PRIME 

0.000330009 

D-6 

0.001243303 

D-6 

0.000355989 

D-6 -0.0008261 

D-6 

CONST 


00000 

NORMAL ECUATICNS 


7.6519783E OC 
1.9707211E 01 
7.7388299E-01 
-1.2076450E 01 


2.3356000E-01 

7.5781000E-02 

-2.97L940CE-01 

-3.729700CE-01 


1.0916790E 00 
-4.6823550E 00 
5.6323100E-C1 
-1.7265100E 00 


3.7076520E OO 
1.0352294E 01 
2. 2502300E-01 
-5.93591S9E 00 


5.3000899E-10 
1.2453030E-09 
3.5598900E-10 
-8.2609999E-1 0 


21 

22 

23 

24 

25 

26 

0.00012539 

-0.00081334 

0.00006408 

-0.00008347 

-0.00000000 

-0.00000000 
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Table C-3 

Listing of CHVL, Subroutines, and a Sample Case. 

$JOB 1203P003 126 M32RK. M3081C1234 
SEXECUTE IB JOB 

SIBJOB SOURCE. MAP »G0 

SIBFTC CHVL DECK.M94.XR7 

C 

DOUBLE PRECISION SA.SE .S02 .SO.S 12 . S I .C02 .CO.SM . SGZ2 , SGZ . P I » TWOPI » 

X RAD » EX »EP2 » EP *SB» T IME.SG.E .SINE* COSE .ESE.ECE.BK* RVZ.RNUM. 

X DELRHO.P.Q.R.SN2.SN,DSQRT.DCOS.DSIN.DELTAE,DELTAA.DELTAG. 

X DELTAI .DELCO.DELSO.DELTAN.SANEW .SENEW .SONEW . SINEW .CONEW . 

X SGZNEW.SNNEW .A.B.C.DELTAT .TZERO.X JD .DTCHEB . SMI NNR » EEK * TAU* 

X EX.DELl.DEL2.DtL3.ASTAR.BSTAR.GGTAR.DlFF 
DIMENSION XG( 3 .1001 J 
DIMENSION A ( 3 ) »B(3) * C ( 3 ) 

DIMENSION H ( 6 ) . W ( 3 ) .RV (3 ) .RNUM ( 3 * 501 ) .D IFF ( 3 .201 ♦ 10 ) 

DIMENSION IT1TLEU2) 

DIMENSION AA( 7 ) .BB( 7 .7 ) .BM( 4 .4 ) .CM ( 4 . 1 ) 

DIMENSION P<3) .0(3) .R<3) . DELRHO ( 3 ) .RVZ ( 3 ) 

1 FORMAT (4D18.1) 

2 FORMAT (26H0INITIAL CONDITION CHANGES/ 1P6E 15. 7 ) 

3 FORMAT (IX) 

4 FORMAT ( 10X3D20 • 10 ) 

5 FORMAT ( I8.6P6F12.4) 

6 FORMAT (3F14.8) 

7 FORMAT < 12HlDIFFERENCtS ) 

8 FORMAT (1PD19.10.10E11.2I 

9 FORMAT (/2(7X1HX1X11HALPHA*1.E102X10HBETa*1.E101X11HGAMMA*1.E10>// 
91) 

10 FORMAT <2( I8.6P3F12.4) ) 

14 FORMAT ( 1H015X5HALPHA7X5HBETA 7X5HGAMMA6x6HALPHA*6X6HBETA* 6X6HGAM 
141MA*7X2HDA10X2HDB10X2HDC ) 

15 FORMAT (1415) 

17 FORMAT (8HK.DELTA TF12.7/6H TZEROF14.7/3H EPF17.7/8H JUL DAYF10.1. 
1718H DTCHEB F9.0) 

18 FORMAT ( F 10. 1 .6P9F12.4 ) 

19 FORMAT ( 2 1H0 JD T AU .SG2 .E25X5HRZER09X1HW 1 1X1HR13X4HRNUM11X5H 

191A.B.C6X8HA*.B*.C*4X8HDA»DB.DC/0PF10.1 »F 12*8 »3F 12. 7 » 1PD18 . 10 .6P3F12 
192.4) 

20 FORMAT ( 10X0PF 12 . 4 . 3F12 . 7 . 1PD18. 10 .6P3F12 .4 ) 

21 FORMAT ( 32HKNEW ELEMENTS FOR DISTURBED BODY) 

23 FORMAT ( 1HL13X2HZ1 14X2HZ214X2HZ314X2HZ414X2HZ5 14X2HZ6/4X6P6F16. 8 ) 

24 FORMAT (22HL CHANGE IN ELEMENTS/ /I 1X2 HSA16X2HSE16X 
X 2HS016X2HSI 16X2HC016X3HSGZ15X2HSN/7E18.8 ) 

27 FORMAT (24H0R00T MEAN SUUARE ERRORS/1H08X5HALPHA10X5H BETA10X5HGAM 
271MA/6P6F15.6) 

28 FORMAT (6P6F12.6) 

29 FORMAT ( 17HONORMAL EQUAT IONS/1X ) 

30 FORMAT (1P7E15.7) 

31 FORMAT ( 24H0MAX IMUM ABSOLUTE ERRORS/ 1H08X5HALPHA10X5H BETA10X5HGAM 
31 1MA/6P6F15 .6 ) 

33 FORMAT (4D18.11) 

36 FORMAT (26H0NTIME IPRINT IELE ID IPRM) 

41 FORMAT (12A6) 

43 FORMAT (1H1) 

44 FORMAT (8HJSMALL A/D24.16) 

45 FORMAT ( 8HJSMALL E/D24.16) 

46 FORMAT ( 12HJSMALL 0MEGA/D24. 16 t 

47 FORMAT (8HJSMALL I/D24.16) 

48 FORMAT (10HJCAP OMEGA/D24.16 ) 

49 FORMAT ( 8HJ SMALL M/D24.16) 

50 FORMAT ( 7HJG ZERO/D24.16) 

51 FORMAT (8HJSMALL N/D24.16) 
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Table C-3 (Continued) 


52 FORMAT ( 1HJ10X1HP19X 1HQ19X1HR 19X1HA19X1HB19X1HC/ ( 1P6D20. 12 ) I 

53 FORMAT ( 7HJSM1 NNR/D24. 16 ) 

PI = 3. 14 1592653589793 

TWOPI=2.*PI 

RAD=TWOPI/360. 

55 CONTINUE 
DO 59 1*1.3 
DO 59 J=1»1Q01 
59 XG(I,J)“0. 

READ (5.41) ITITLE 

READ (5,15) NTIME, IPRINT, IELE, ID, IPRM 
C 

C NTIME IS NUMBER OF TABULATED DATES 
C IPRINT CONTROLS PRINT 

C I ELE = 0 FOR REGULAR RUN(2 ITERATIONS. . ■ NO. OF ITERATIONS FOR OTHER RUNS 
C ID GREATER THAN ZERO FOR TEST OF COORDINATE DIFFERENCES 
C IPRM = 0 FOR SUN AS PRIMARY. GREATER THAN ZERO FOR EARTH PRIMARY 

C 

READ (5,33) SA »SE »S02 .SI2.C02.SM.SGZ2 »SM j NNR 
READ (5.1) DELTAT » TZERO »XJD » DTCHEB 
C TZERO IS STARTING TIME FOR EVALUATION COUNTED IN DAYS FROM EPOCH 
C DELTAT IS STEP INTERVAL FOR EVALUATION 
C XJD IS JULIAN DATE OF TZERO 

C DTCHEB IS TOTAL INTERVAL OF VALIDITY OF SERIES IN DAYS 
READ (5.28) Z1 ,Z2 ,Z3 ,Z4 ,Z5 ,Z6 
S0=S02*RAD 
SI =SI2*RAD 
C0=C02*RAD 
SGZ=SGZ2*RAD 
EK=0 .017202098950+00 
EP2=23.4457888889 
EP=EP2*RAD 

C IF ( 1PRMJSUN, SUN, EARTH 
IF (IPRM) 72,72,70 
70 EK=4 .4619969180+0 
EP'0.0 
72 CONTINUE 

SN=EK*DSQRT ( 1 . +SM 1 NNR+SM ) /SA**1 .5 
SN2=SN/RAD 
WRITE (6,43) 

WRITE I 6 ,41 ) I T 1 TLE 
WRITE (6,36) 

WRITE (6.15) NTIME. IPRINT, IELE, ID, IPRM 
WRITE (6.44JSA 
WRITE (6.45ISE 
WRITE (6,46) S02 
WRITE ( 6 , 47 ) S I 2 
WRITE (6,48)002 
WRITE (6,49 ) SM 
WRITE (6.50ISGZ2 
WRITE (6.51ISN2 
WRITE (6.53) SMI NNR 

WRITE ( 6, 17 IDELTAT, TZERO, EP2, XJD. DTCHEB 
WRITE (6»23)Z1»Z2»Z3»Z4.Z5,Z6 
WRITE (6,43) 

WRITE (6.41) ITITLE 

READ CHEbYCHEV SERIES PERTURBATION TERMS 


WRITE (6,9) 
NTRM-0 
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Table C-3 (Continued) 


90 READ <5. 5) K. ( H( I ) . 1=1 ,6 ) 

KP=K+1 

WRITE (6.10) K.*(H<II*I = 1*3)*K.P»(H(I)»I=4*6) 
IF (K-1001) 100.120*120 
100 IC1 = K+1 
KP=K+2 
NTRM=KP 
DO 110 1=1*3 
J*l + 3 

XG ( I *K1 ) =H ( I ) 

XG(I.KP)=H(J) 

110 CONTINUE 
GO TO 90 
120 CONTINUE 

READ GENERAL PERTURBATIONS TERMS 
READ COORDINATES 


CARDS REMOVED TO MODIFY FOR RES. 3 BOY. PROB* 

160 CONTINUE 
165 CONTINUE 
I TRMAX=2 

IF ( IELE ) 184,184,182 
182 ITRMAX= IELE 

184 CONTINUE 

ENTRY FOR COMPARISON 

ITRATN=0 

185 CONTINUE 
SB=SA*DSQRT ( l.-SE**2 ) 

PORDP2 

CALL PORDP2 (SO *SI *CO*P*Q*R* EP ) 

DO 186 1=1*3 
A ( I ) =SA*P ( I > 

B( I )=SB*Q< I ) 

C ( I )=SA*R( I ) 

186 CONTINUE 

WRITE (6.52) (P( I ) »Q( I ) *R( I ) »A( I )»B( I )*C( I )»I=1*3) 

ASSIGN 322 TO NS1 
188 WRITE (6*43) 

WRITE (6.41) IT1TLE 
ITRATN=ITRATN+1 
IF (IPRINT-D192. 192*190 
190 ASSIGN 323 TO NS1 
GO TO 194 
192 WRITE (6.14) 

194 CONTINUE 

ZERO OUT BB(I.J). 

DO 196 1=1.7 
DO 196 J*l*7 
196 BB( 1 *J)=0. 

AAA=0. 

BBB=0. 

GGG=0. 

DAMX=0. 

DBMX=0. 

DGMX=0. 
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Table C-3 (Continued) 


EQUALLY SPACED TIME INTERVALS DO LOOP. 


NTR2=NTRM+1 
NTR =NTRM-1 
EEX=2./DTCHEb 
DAY=XJU-DELTAT 
T IME=-DELTAT+TZERO 
DO 324 N=1 *NT IME 
T IME=T IME+DELTAT 
TAU=EEK*T IME 
DAY=DAY+DELTAT 
SNT =SN*T IME 
SG=SGZ+SN*T IME 
SG2=SG/RAD 

COMPUTE ECCENTRIC ANOMALY. 


E=SG 

SINE=DSIN(SG) 

COSE=DCOS(SG) 

SING=SINE 

COSG=COSE 

DO 198 M= 1 » 3 

ESE=SE*S I NE 

ECE=l.-SE*COSE 

BK = ( SG-E+ESE ) /ECE 

E=E+BK-.5*bK*BK*ESE/ECE 

COSE=DCOS ( E ) 

S I NE=DS IN ( E ) 

198 CONTINUE 

COS2E=l.-2.*SINE**2 

SIN2E=2.*SINE*COSE 

Cl=l./(l.-SE*COSE) 

E2=E/RAD 

ALPHA=0. 

BETA=0. 

GAMMA=0. 

EVALUATE CHEBYCHEV SERIES 
EX=2.*TAU 

IF (NTR) 240.200*210 
200 ALPHA=.5*XG( 1.1 ) 

BETA = . 5*XG (2.1) 

GAMMA= . 5*XG (3.1) 

GO TO 240 
210 CONTINUE 

DO 230 1=1.3 
B2=0. 

81 = 0 • 

DO 220 K=1»NTR 
I1=NTR2-K. 

bZ=EX*bl-b2+XG( I . 1 1 ) 

B2=bl 

B1=BZ 

220 CONTINUE 

H ( I )=.5*(EX*B1-2.*B2+XG( 1.1) ) 
230 CONTINUE 
ALPHA=H ( 1 ) 

BETA =H ( 2 ) 
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Table C-3 (Continued) 


6AMMA=H < 3 ) 

240 CONTINUE 

EVALUATE GENERAL PERTURBATIONS SERIES 


ADO IN CONSTANTS OF INTERGRAT ION. 

ALPHA=ALPHA+2 . *Z3+Z1* ( COSE-SE ) +Z2*S INE 

BETA=BETA+Z4+Z1*( 3.*SE*SNT- ( 2 .-SE**2 ) *S I NE+. 5*SE*SIN2E » + 

X Z2*(2.*COSE-.5*SE*COS2E)+Z3*(-3.*SNt) 

GAMMA=GAMMA+Z5* < COSE-SE ) +Z6*S I NE 

COMPUTt PERTURBATION VECTOR. 

DO 295 1=1.3 

RVZ( I )=PI I ) *SA* ( COSE-SE ) +Q ( I )#SB*SINE 
W(I)=C1*(**SINE*A( I )+COSE*B< I ) ! 

RV ( I ) = ( 1 • +ALPHA ) *RVZ ( I )+BETA*W< I ) +GAMMA*c I I ) 

295 CONTINUE 

CARDS HAVE REMOVED AND ADDED TO MODIFY FOR RES. 3 BDY. PROB 
ASTAR=4081 .5946 E-6 
BSTAR=Q. 

GSTAR=0. 

COMPUTE BB(I.J). 

AA ( 1 ) =ALPHA-AST AR 
AMA=AA ( 1 ) 

DAMX = AMAX 1 ( DAMX • ABS ! AMA ) ) 

AAA=AAA+AMA*AMA 
AA ( 2 ) =COSE-SE 
AA ( 3 I =S INE 
AA(4)=2. 

DO 310 J= 1 *4 
DO 310 I=J.4 

310 BB( I * J ) =BB ( l.J) +AA ( I ) *AA ( J ) 

AA( 1 )=BETA“BSTAR 
BMB=AA(1) 

DbMX=AMAXl (DBMX.ABSI BMb ) ) 

Bbb=BBb+bMB*bMB 

AA(2)=3.*SE*SNT-(2.-SE**2 )*SINE+.5*SE*SIn2E 
AA ( 3 ) =2 .*COSE-.5*SE*COS2E 
AA(4)=-3.*SNT 
AA( 5 ) = 1 . 

DO 315 J=1 ♦ 5 
DO 315 I = J . 5 

315 BB( I * J ) =Bb ( I » J >+AA( I I *AA ( J ( 

AA( 1 ) =GAMMA-GST AR 
GMG=AA ( 1 ) 

DGMX=AMAX1(DGMX*ABS(GMG) ) 

GGG=GGG+GMG*GMG 
AA( 6 ) =COSE-SE 
AA ( 7 ) =S INE 

BB( 1 .1 ) =BB( 1.1) +AA C 1 ) *AA ( 1 ) 

BBI6.1 )=BB(6.1 )+AA(6)*AA( 1 ) 

BB ( 7 . 1 ) =BB (7.1) +AAI 7 ) *AA ( 1 ) 

DO 320 J=6 . 7 
DO 320 I = J . 7 

320 BB ( I . J ) =BB( l.J) +AA! I )*AA ( J ) 

ATA=ALPHA 


55 



n> n n 


Table C-3 (Continued) 


BTA= BETA 
GTA=6AMMA 
ATS=ASTAR 
BTS=BSTAR 
GTS=GSTAR 
BAT = AMA 
BBT = BMB 
OGT= GMG 

GO TO NS1.( 322.323) 

322 WRITE (6,18) DAY » ATA,BTA,GTA,ATS,BTS,GTS,DAT ,DBT ,BGT 
GO TO 324 

323 CONTINUE 

WRITE (6,19) DAY»TAU.RVZ(1),W(1),RV(1),RNUM(1,N), ATA.ATS.DAT 
WRITE (6.20) SG2,RVZ(2),W(2),RV.(2),RNUM(2,N),BTA,BTS,DBT 

WRITE (6.20) E2.RVZ(3),W(3),RV(3),RNUM(3,N),GTA,GTS,DGT 

324 CONTINUE 

DO 325 J=3 « 5 
IJ=J-1 

DO 325 1=2. IJ 

325 Bb< I , J ) =BB( J , I ) 

B6 ( 6 , 7 ) =BB (7,6) 

DO 330 1=1.4 

CM ( 1,1) =— BB ( 1+1,1) 

DO 330 J= 1.4 
330 BM( 1 « J ) =BB( 1+1 » J+l ) 

WRITE (6,43) 

WRITE (6.41) ITITLE 
WRITE (6,29) 

WRITE (6,30) ( (BB( I , J ) , I =2 ,7 ) ,BB ( J , 1 ) » J=2 , 7 ) 

CALL MATINV(BM»4, CM, 1, DETERM) 


Z1=CM( 1,1) +Z1 
Z2=CM (2,1) +Z2 
Z3=CM (3,1) +Z3 
Z4=CM (4,1) +Z4 


Z5= ( BB l 7 , 7 ) *BB (6,1) -Bb (7.1) *BB (6,7) ) / < BB (7,6) *BB (6,7) 

X -BB(7,7)*BB(6.6) ) +Z5 

Z6= ( BB ( 7 , 1 ) *bB (6,6) -BB ( 6 . 1 > *bB ( 7 » 6 ) ) / ( BB ( 6 » 7 ) *B6 < 7 ,6 ) 

X — BB ( 7 , 7 ) *BB (6,6) ) +Z6 

WRITE (6,23) Z1,Z2,Z3,Z4,Z5,Z6 


COMPUTE CHANGE IN ELEMENTS. 


C0N1=SQRT ( 1 .-SE**2 ) 

DELTAP=Z6/C0N1 
DELTAQ=-Z5 
DELTAR=-C0N1*Z2/SE 
DELTAE=-( l.-SE**2 )*Z1 
DELTAA=SA*(2.*Z3-2.*SE*Z1 ) 

DELT AG=Z4+ ( 2 . +SE**2 ) *Z2 / ( 2 . *SE ) 
COSSO=COS(SO) 

S I NSO=S I N ( SO ) 

DELTAI =DELTAP*COSSO— DELTAO*S INSO 
COSSI=COS(SI ) 

SINSI=SIN(SI ) 

DELCO= ( DELT AP*S I NSO+DEL T AQ*COSSO ) / S I NS I 

DELSO=DELT AR— COSS I *DELCO 

SANEW=SA+DELTAA 

SENEW=SE+DELTAE 

SONEW=SO+DELSO 

S INEW=SI+DELTA I 

CONEW=CO+DELCO 
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Table C-3 (Continued) 


SGZNEW=SGZ+DELTAG 

SNNEW=ER*DSQRT ( 1 . + SM I NNR+SM i / SANEW** 1 . 5 

DELTAN=SNNEW-SN 

DELTAI=DELTAI/RAD 

U£LSO=UELSO/RAD 

DELCO=DELCO/RAD 

DELTAG=DELTAG/RAD 

DELTAN=DELTAN/RAD 

S 1 2 = S I NEW /RAD 

SU2=SONEW/RAD 

C02=C0N£W/RAD 

SGZ2=SGZNEW/RAD 

SN2=SNNEW/RAD 

ENT IME=NT 1ME 

AAA=SQRT ( AAA/ENT IMt ) 

BBB=SQRT ( BBB/ENT 1ME ) 

6GG=SQRT ( GG6/ENT IME ) 

WRITE (6*31 > DAMX . DbMX . UGMX 

WRITE (6*27 ) AAA * BBd *GGG 

WRITE ( 6 * 24 ) UELTAA.DELTaE.DELSO.DELTAI . 

X DELCO.DELTAG.DtLTAN 


WRITE 

(6.21 ) 

WRITE 

(6.44) SANEW 

WRITE 

(6.45 JSENEW 

WRITE 

(6.46JS02 

WRITE 

(6.471SI2 

WRITE 

(6.481C02 

WRITE 

(6.49JSM 

WRITE 

(6.501SGZ2 

WRITE 

<6.51 )SN2 

WRITE 

(6.53ISMINNR 


COMPUTE CHANGES IN INITIAL CONDITIONS 

SG=SGZ 
E=SG 

S I NE=DS I N ( SG ) 

COSE=DCOS ( SG ) 

DO 340 M=1 *3 
ESE=SE*S I NE 
ECE=1.-SE*C0SE 
6K.= ( SG-E+ESE ) /ECE 
E = E+BK-.5*BK*8K.*ESE/ECE 
C0SE=DC0S(E) 

SINE=DSIN ( E ) 

340 CONTINUE 

C0S2E=1 .-2.*SINE**2 
SIN2E=2.*SINE*COSE 
DA=2.*Z3+Z1*( COSE-SE ) +Z2*SINE 
DG= Z5* ( COSE-SE ) +Z6*SINE 

DB=Z4+Zl*(-(2.-SE*SE)*SINE+.5*SE*SIN2E)+z2*(2.*COSE-.5*SE*COS2E) 
ADR=1./ ( l.-SE*C0SE) 

DAP=ADR*(-Z1*SINE+Z2*C0SE ) 

DGP=ADR* ( -Z5*S INE+Z6*C0SE ) 

DBP=3. * ( St*Z 1-Z3 ) +ADR* ( ( - ( 2 . -SE*SE ) *C0SE+SE*C0S2E ) *Z 1 
1 +(- 2. *SINE+SE*SIN2E)*Z2) 

WRITE (6*2) DA *Db »DG * DAP »DBP » DGP 
IF ( ITRATN-ITRMAX ) 188.55.55 
END 

SIBFTC PQRDP2 DECK.M94.XR7 

SUBROUTINE P0RDP2(S0,SI *CO*P *Q*R*EP ) 
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Table C-3 (Continued) 


C COMPUTE VECTORS P.U.R 

UOUbLE PREc 1 S I ON SO.SI.CO.P.Q.R.EP.EPC.EpS.CCO.SCO.CSO.SSO. 
X CSI .SSI .TP1.TP2.DCOS.DSIN 

DIMENSION P(3) .0(3) * R ( 3 ) 

C EP IS OBLIUUITY OF ECLIPTIC 
EPC = DCOS( EP ) 

EPS=DS I N ( EP ) 

CCO=DCOS ( CO ) 

SCO=DS I N ( CO ) 

CSO=DCOS ( SO ) 

SSO=DS I N ( SO ) 

CSI=DCOS< SI ) 

SSI =D3 I N ( S I ) 

P( 1 )=-CSI*SSO*SCO+CSO*CCO 

TP1 =+CSI*SSO*CCO+CSO*SCO 

TP2 =+SSl*SSO 

P(2 )=EPC*TP1-EPS*TP2 

P(3 )=EPS*TP1+EPC*TP2 

Q(1 )=-CSI*CSO*SCO-SSO*CCO 

TP1 = CSI*CSO*CCO-SSO*SCO 

TP2 = 5S1*CS0 

Q(2 ) =EPC*TP1-EPS*TP2 

Q(3)=EPS*TP1+EPC*TP2 

R ( 1 ) = SSI*SCO 

TP1 =-SSI*CCO 

TP2 = CSI 

R12 )=EPC*TP1-EPS*TP2 

R(3)=EPS*TP1+EPC*TP2 

RETURN 

END 

SIbFTC MATINV DECK. . M94 * XR 7 

SUBROUTINE MAT 1 NV ( A . N . B *M .DET ERM ) 

C 

DIMENSION IPIVoT (4) » A ( 4 *4 ) » B ( 4 * 1 ) . INDEX (4*2 ) *P I VOT ( A ) 

INI T IALIZAT1UN 

10 DETERM=1 .0 
15 DO 20 J=1.N 
20 I P I VOT ( J ) =0 
3U DO 550 I = 1 » N 

SEARCH FOR PIVOT ELEMENT 

AO AMAX=0 • 0 
45 DO 105 J=1*N 
50 IF ( IPIVOT(J)-l) 60*105.60 
60 DO 100 K=1*N 

70 IF ( IP1VOTIK1-1 ) 80.100*740 
80 IF (ABS (AMAX)-ABS (A(J.K))) 85.100.100 
85 I R0W= J 
90 I COLUM=K 
95 AMAX=A(J.K) 

100 CONTINUE 
105 CONTINUE 

110 IPIVOT 1 I COLUM ) = I P I VOT ( ICOLUMJ + l 

INTERCHANGL ROWS TO PUT PIVOT ELEMENT ON DIAGONAL 

130 IF ( IROW-ICOLUM) 140.260.140 
140 DETERM=-DETERM 
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Table C-3 (Continued) 


150 DO 200 L=1 *N 
160 SWAP=A ( I ROW »L ) 

170 A ( I ROW * L ) =A i 1C0LUM » L ) 

200 A ( I COLUM* L ) =SWAP 
205 I F ( M ) 260 ,260, 2 10 
210 DO 250 L=1»M 
220 SWAP=B ( 1 ROM » L ) 

230 B ( I ROW » L ) *B ( i COLUM » L ) 

250 B ( I COLUM » L } = SWAP 

260 INDEX (1*1) = IROW 

270 INDEX! I * 2 ) “ICOLUM 

310 PIVOT ( I )=A( ICOLUM* ICOLUM) 

DIVIDE PIVOT ROW BY PIVOT ELEMENT 

330 A( ICOLUM, IC0LUM)=1.0 
340 DO 350 L=1»N 

350 A ( 1 COLUM* L ) *A( ICOLUM, L) /PIVOT ! I ) 

355 I F ( M ) 380,380,360 
360 DO 370 L= 1 »M 

370 B < ICOLUM, L ) =B! I COLUM, L) /PI VOT ( I ) 

REDUCE NON-PIVOT ROWS 

380 00 550 L1=1,N 

390 IF(Ll-ICOLUM) 400,550,400 

400 T=A( LI, ICOLUM) 

420 A ( L 1 , I COLUM ) =0,0 
430 DO 450 L= 1 »N 

450 A(L1,L)=A(L1«L) -A ( ICOLUM,L)*T 
455 IF ( M ) 550 , 55U , 460 
460 DO 500 L=1»M 

500 B(L1,L)=B(L1,L)-B( ICOLUM,L)*T 
550 CONTINUE 

INTERCHANGE COLUMNS 

600 DO 710 1=1, N 
610 L=N+1— I 

620 IF ( INDEX(L,1)-INDEX(L,2) ) 630,710,630 
630 JROW= I NDEX ( L, 1 ) 

640 JCOLUM= INDEX! L ,2 ) 

650 DO 705 ,K=1»N 
660 SWAP=A ! K, JKOW ) 

670 AIK , JROW J = A (X , JCOLUM ) 

700 A (K, JCOLUM )=SWAP 
705 CONTINUE 
710 CONTINUE 

DO 800 1= 1 »N 
J=N+1 - I 

800 DETERM=DETERM*PIVOT( J) 

740 RETURN 
ENO 
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